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CHAPTER I 


INTRODUCTION 

The atmosphere is an inhomogeneous non-isotropic media which is 
usually in a state of turbulence. The atmosphere is characterized by 
its temperature, wind velocity, and humidity. The variations of these 
parameters comprise non-stationary random processes and as a result the 
fluctuation of the index of refraction of the media is also a non- 
stationary random process. Electromagnetic radiation at optical wave 
lengths propagating through the atmosphere will be greatly affected by 
the random fluctuation in the index of refraction. This causes a very 
complicated scattering to occur which results in amplitude and phase 
variations of the wave. Clearly, these fluctuations are also of a random 
nature . 

The distortion induced by the atmosphere on the propagating wave is 
of great concern in the development of optical tracking and communication? 
systems since performance 'can .be seriously^/degraded due to this effect. 
For example, amplitude fluctuations cause the signal- to-noise ratio to 
be reduced in incoherent detection systems. Atmosphere distortion tends 
to reduce the coherence of the wave which in turn reduces the effective 
power level of the signal. This will also decrease the signal-to— noise 
ratio. Loss of coherence is a problem in systems utilizing coherent de- 
tection where heterodyne action must be achieved. Wave front tilt due 
to phase variations of the wave induces errors in optical tracking re- 
ceivers since they view this tilt as an apparent angle of arrival. 

An important aspect in the design of optical systems that is de- 
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pendent on a knowledge of atmosnheric distortion is the size of the optical 
components. Atmospheric effects can be somewhat • overcome by employing 
large receiving optics which will tend to average out the wave fluctuations. 
This advantage is limited by the high expense and difficulty in fabrication 
of such components. The designer must then choose optimum size components 
which require that he have a thorough knowledge of the atmospheric pro- 
blem. It is clear that there is a pressing need for an accurate 
mathematical model of the atmosphere. 

Statistical methods must be employed to analyze this problem since 
the wave fluctuations are random processes. The first step in developing 
a statistical model is to determine the probability density function of 
the random process. Theoretical considerations have predicted a log- 
normal distribution for the amplitude and normal for the phase. This 
theory has been experimentally verified for visible wave lengths, but re- 
sults of current investigations in the infrared region of 10.6 microns 
have been inconsistent. 

D. L. Fried^ has made scintillation measurements at this wave length 
over a 1 km. range using a point detector. His results do not confirm 
the hypothesis that intensity scintillation is log-normally distributed. 

He suggests that this may be a genuine feature of 10.6 micron scintilla- 
tion but draws no definite conclusion since detector noise and nonlinearity 

problems in taking measurements could have influenced his results . 

2 

Richard Kerr of the Oregon Graduate Research Center has conducted 
multiwave length laser propagation studies over a mile path and claims 

O 

confirmation of log-normal statistics for wave lengths of 4880 A and 10.6 

3 

microns. In addition Fitzmaurice, Bufton, and Minott have also con- 
cluded that scintillation at 10. 6; micron fits the log-normal model. Their 
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work was done over a 2.4 km. path. Both investigators used point 
source detection. 

The effect of the atmosphere on 10.6 micron propagation is important 
since the popular laser emits radiation at this wave length. This 
type laser is attractive for application in. optical systems due to its 
high efficiency and high output power capability. In addition, the 
atmospheric effect at this wave length is much less than that at visible 
wave lengths . 

The purpose of this study is experimentally to investigate the 
statisitical properties of scintillation and the signal-to-noise ratio 
of heterodyne detection for a CO^ laser beam propagated over a 3.2 km. 
path. Both scintillation and heterodyne measurements have been made for 
a variety of receiving 'aperture sizes ranging from two to ten cm. 

A brief discussion of the theory which is referred to in current 
literature is presented in Chapter II, The necessary statistical con- 
cepts are introduced before a qualitative description of atmospheric 
turbulence is ^iven. Finally, the physical significance of aperture 

i 

averaging is discussed. 

Chapter III gives a detailed description of the experiment. De- 
scribed is the equipment, its alignment and check out as well as a 
discussion ore the techniques used 'to make the measurements. 

The handling and reducing of the data is given in Chapter IV. This 
includes a discussion on the conversion of analog data to digital form 
for direct .use on a digital computer. An outline of the computer pro- 
gram which reduces the data is presented. The theory used to calculate 
> 

aperture effects is also given. 

Chapter V is concerned with interpreting the reduced data to de- 
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termine if the hypothesis of the log-normal distribution for intensity 
scintillation is valid for this wave length. This chapter also in- 
cludes results of calculations for the refractive index structure con- 
stant with and without aperture averaging corrections. 

Chapter VI contains the summary and conclusions of this study as 
well as recommendations for further study. A complete documentation 
for the computer program is given in the' Appendix. 



CHAPTER II 


THEORY 


A. Statistical Concepts 

It is necessary to give a discussion on pertinent statistical con- 
cepts as a prelude to presenting .a- qualitative discussion on the 
theoretical aspects of- the atmospheric problem. 

The random processes are described in terms of parameters which are 
random variables. The value of any such function at- a fixed instant of 

time is a random variable having definite probability density- function. 

1 

The process may further be described by its auto-covariance function at 
times t^ and t 2 

AC{f( tl ), f(t 2 )} = - <^(tj> ][f(t 2 ) - <f(t 2 )>]> 2-1 

\ 

where indicates an ensemble average. The auto-covariance function 
reduces to' the correlation function 

b[f (t-^) 9 f(t 2 >] = <^(t 1 )f(t 2 ^> 2-2 

for processes where the mean value is zero. The auto-covariance- function 
characterizes the mutual relation between the fluctuations at different 
instants of time. The mean value of the random variable can be a con- 
stant or can change Ttfith time. Similarly, the auto-covariance function 
can either depend only on the difference between the times t^ and t 2 or 
else it can depend on the positions of the points on the time axis. The 
first case would occur when the statistical relation between the 
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fluctuations of the variable at different instants of time does not 
change with time. A random function is called stationary if its mean 
value does not depend on time and its auto-covariance function depends 
only on the difference between observation times. 

The mean value of the meteorological parameters of the "atmosphere 
such as temperature, wind velocity, and humidity undergo comparatively 
slow and smooth changes. These variables are non-stationary processes 
if the definition of stationarity is strictly applied. It is difficult 
to determine which changes in the fluctuation are to be regarded as slow 
changes in the mean and which are to be regarded as slow fluctuations of 
the function. 

To avoid this difficulty and to describe random functions which 

have the above characteristics, the structure function is used instead 

of the correlation function. This function was first introduced by 

4 5 

Kolmogorov * . The basic idea behind this method is to use the difference 
function 


F T (t) = f(t+r) - f(t) 2-3 

instead of the non-stationary function f(t). For values of x which are 
not too large, slow changes in the function f(t) do not affect the 
value of the difference function which means that it can be considered a 
stationary • random function. The function f(t) is called a random function 
with stationary increments. To derive an expression for the structure 
function consider the transformation of the correlation function for 
F(t ) and F(t„) : 

T 1 T ^ 


B(t i’ V - 


2-4 



B (t-L> t 2 ) = <ff(t 1+ t) - f (t x )] [f (t 2 +r) - f(t 2 )£> 
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Using the algebraic identity 

(a-b) (c-d) = j [ (a-d)^ + (b-c) 2 - (a-c) 2 - (b-d) 2 ] 2-6 


we have 






[f(t 1 +x) - 


f<t 2 )3 






{f<t 1 +r) - f<t 2 -H)]‘ 


1 

2 


[f.(t 1 ) - f<t 2 )r 
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Thus the correlation function is expressed as a linear combination of 
functions of the form 


W 



< 


[f(t ) ^ f(t )]' 
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which is called the structure function of the random process. The form, 
of the structure function more commonly used 


V T) 


[f(t+T> - £{t)] 



2-9 


is the basic characteristic of a random process with stationary in- 
crements. The value of D(t) characterizes the intensity of those 
fluctuations of f(t) which are smaller than or are comparable with t. 


B. Nature of Atmospheric Turbulence 

The statistical theory of turbulence was initiated in the works 
.of Friedmann and Keller^. This theory was greatly amended in 1941 

g 

when A. N. Kolmogorov and A. M. Obukhov established the laws which 
characterize the basic properties of the microstructure of turbulent 
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flow at very large Reynolds numbers. The following discussion is drawn 

g 

from Tatarski's work on the Kolmogorov theory. 

Consider the atmosphere to be a viscous fluid in a state of 
laminar flow. This flow can be characterized by its viscosity v, velocity 
v, and the characteristic length L. The quantity L characterizes the 
dimensions of the flow as a whole and arises from the boundary con- 
ditions of the fluid dynamics problem. This laminar flow will be stable 
if the Reynolds number 



does not exceed a certain critical value. 

Suppose that for some reason a velocity fluctuation occurs in a 
region of size Z of the initial laminar flow. The value of Reynolds 
number will increase and the laminar motion will lose stability. The 
result of this instability is the formation of a secondary flow or eddies 
within L which will have their own Reynolds number R . As the Reynolds 
number for the overall flow is increased, R^ will increase causing the 
secondary flow to break up into smaller scale eddies. These new eddies 
now give energy to even smaller eddies and the process continues until 
an eddy with a Reynolds number less than the critical value is formed. 

The atmospheric turbulence can be considered as consisting of many 
circulating eddies having different flow characteristics. Eddies are 
usually described in terms of an inner and outer scale of turbulence. 

These are measures of the characteristic sizes of the smallest and largest 
eddies which exist at the time of interest. Figure 1 may aid in visualizing 
this process. The outer scale is the physical dimension of the largest 
eddy. The inner scale of turbulence is roughly the size of the smallest 




Figure 1. Visualization of Microstructure of Turbulence. 
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stable eddy. It can be more precisely defined in terms of a characteristic 
of the longitudinal velocity structure function 


D 

rr 


<‘ 


CV< ri ) - V(r 2 )] 
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where V(r^) is the projection of the velocity at the point r^ along the 
direction of r, and V(r 2 ) is the same quantity at the point 


For r << a 

o 


and for r >> £ 
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r l + r 
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2-13 


For r on the order of L the structure function saturates. The inner 

scale of turbulence is then defined mathematically as the value for 

D where the functions in equations 2-12 and 2-13 intersect, 
rr 

Each eddy or cell in the field of turbulence can be considered 
locally isotropic and homogeneous, and as a result it will have a certain 
index of refraction, which we will assume to be constant throughout the 
cell. The index of refraction will in general differ from cell to cell. 
As an electromagnetic wave passes through each cell two things occur: 
First, the phase of the wave is advanced or retarded in a random manner 
due to the index of refraction of the cell. Secondly, the wave is 
scattered due to the interfaces between the cells. This causes the in- 
tensity of the beam to be distributed at random after traveling through 
a large number of cells . 

From the model of the atmosphere just developed some insight as to 
the statistical characteristics of phase and amplitude variation can be 
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gained. Since the random variations in phase add to each other, the 
Central Limit Theorem^ can be applied to predict that the distribution 
of the phase across the wave front is normal. The amplitude of the wave 
has a distinctly different character. ^After each refraction the intensity 
is the product of the intensity before scattering and the variation due 
to refraction. The intensity variations are then due to a product of 
probabilities. If the Central Limit Theorem is applied to the logarithm 
of the intensities, they will be normally distributed. This leads to a 
log-normal probability density function for the amplitude distribution. 
Because of this, it is customary to describe waves propagating through 
the atmosphere in terms of their phase <j>(r,t) and log-amplitude L (r,t), 

Ql 

where L (r,t) is given by 
3 . 


L a Cr ' t) * 
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or in terms on intensities 




2-15 


where A and I are mean values. Using these equations the complex repre- 
sentations of the wave becomes 


A(r,t) exp [L(r ,t) + j<j>(r,t)] 


■ 2-16 


C. Aperture Averaging 

Collection of light with large aperture optical systems tends to 
average out atmospherically induced intensity and phase fluctuations . 
This causes a smaller variance for both and alters the statistical prop- 
-erties of the intensity variation. 
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Before discussing the principles of aperture averaging it is 
necessary to define correlation distance r . Consider the intensity at 
two points separated by a distance of r. For r equal zero the correlation 
function will be unity. As r increases, the correlation function decreases 
with zero as its lower bound. How rapidly the correlation function de- 
creases with increasing r, is related to the strength of atmospheric 
turbulence, r^ is defined as the distance at which the intensity variations 

of the two points are no longer correlated, r could also be defined as 

o 

the distance at which the variations become statistically independent. 
r Q usually varies from a few millimeters for very strong turbulence to 
several centimeters for mild turbulence. 

Consider the aperture shown in Figure 2 to be the aperture of an 
optical system which collects and focuses light from a diverging beam. 

Let the light intensity across the illuminated aperture be divided into 
n finite circles of radius r^. The intensity within these circles will 
be assumed to be highly correlated. The collection and focusing process 
can be thought of as adding the intensity contribution of each circle 
on the aperture in the focal plane. Averaging of the intensity fluctua- 
tions of each circle across the aperture will occur at the focus if the 

diameter of the aperture is much larger than r such that it contains 

o 

many circles. 

The intensity fluctuation in the circles of radius r are log- 
normally distributed. Since the aperture adds the intensities, it will 
also add the probability density functions . The Central Limit Theorem 
can then be applied to predict that the intensity variations at the focus 
should be normally distributed. In order to apply the Central Limit 
Theorem, we must assume that the average intensities of the circles are 
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of the same order of magnitude across the aperture and that the variance 
from circle to circle should not change significantly. 



CHAPTER III 


EXPERIMENTAL 


A. Description of Equipment 

The experimental measurements were made during the S umm er of 1969 
at the Marshall Space Flight Center's optical range located on Redstone 
Arsenal near Huntsville, Alabama. 

The transmitting system was located in an astronomical observatory 
on the crest of Madkin Mountain. The receiving and recording systems 
were located in the Astrionics Laboratory complex in a special building 
equipped with a large mirror periscope so that the optical equipment 
could be conveniently placed on the ground level yet have a clear optical 
path to the mountain. The height of the periscope was about 15 feet 
above ground level. The optical path extended in a southwesterly 
direction for a distance of 3.2 km. The transmitter was about 220 meters 
above the receiver so that the optical path was at an angle of 4° with 
the horizontal. Except for a few small buildings and a parking lot 
paved with bituminous material, the optical path lay mostly over wooded 
terrain. 

The transmitter and receiver were constructed by Minneapolis 
Honeywell Corporation for Marshall Space Flight Center and have been 

g 

described in the literature . The transmitter consisted of a 5 watt CO^ 
flowing gas laser with a 10 cm. off axis cassegrainian collimator as 
shown in Figure 3. The laser was designed to have good short and long 
term frequency stability. This was accomplished in part by constructing 
the cavity of the low expansion material cervit, which has an expansion 
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coefficient of less than 1 x 10 ^ (1/°C). To further stabilize the 
laser, water held at a constant temperature to better than 0.1°C was 
continuously circulated through the cervit yoke. The laser in the trans- 
mitter was frequency modulated by applying the modulation voltage to a 
piezoelectric cylinder on which one of the end mirrors was mounted. The 
other end mirror consisted of an Irtran output coupler that was also 
attached to a piezoelectric cylinder, which provided laser transition 
selection, A DC bias was applied to the cylinder to select the desired 
transition. In addition, the transmitter included a mechanical chopper 
that was originally intended to be used for alignment purposes. The 
transmitter unit also contained the necessary electronics to produce a 
modulation voltage for both carrier or direct modulation. A block diagram 
of the transmitting unit is given in Figure 4. 

* 

The receiving unit was housed in a cabinet identical to that of the 
transmitter as shown in Figure 5. The receiver consisted of a 10 cm. 
off axis cassegrainian telescope, a local oscillator laser identical to 
that of the transmitter, combining optics, and a mercury doped cadmium 
telluride detector. This is an alloy detector having a spectral re- 
sponse between 8 and 14 microns. The detector was cryogenic and required 
an operating temperature of 77°K, which was obtained by using liquid 
nitrogen. The operation of the receiving unit can be described with 
reference to Figure 6. The transmitter and the local oscillator signal 
are made spatially colinear by means of a germanium beam splitter and 
combined on the surface of the mercury doped cadmium telluride detector. 
The local oscillator frequency is offset by 10 Mhz from the transmitting 
laser. The 10 Mhz beat frequency produced by the detector is amplified 
with a 10 Mhz center frequency, intermediate frequency amplifier. This 
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amplifier has a 2 Mhz bandwidth and a 110 db gain. The intermediate 
frequency amplifier is followed by a limiter that eliminates amplitude 
variations. A 10 Mhz discriminator provides an error. signal for the 
local laser feedback loop which consists of a low pass filter and a 
feedback amplifier. The feedback amplifier drives a piezoelectric cylinder 
on which one end mirror of the laser is mounted. This provides auto- 
matic frequency control of the local oscillator laser. 

The data acquisition system was located at the receiver terminal 
and consisted of an Ampex 14 channel analog F.M. tape recorder, a variable 
gain AC amplifier with good low frequency response, and two oscilloscopes 
used for monitoring purposes. Figure 7 gives a block diagram of this 
system. A spectrum analyzer was also available to check the 10 Mhz beat 
signal in the receiver. 

The monitoring of both the input to the amplifier and the recorder 
was necessary to insure that they were operating within their linear 
range. Especially critical was the input level to the recorder, since 
its linear range for input voltage was ±1 volt. To be safe we operated 
within ±.5 volt. 

B. System Alignment 

It was necessary to measure the laser and amplifier noise of the 
system to insure that it would not have a significant effect on 
measurements made through the atmosphere. Noise measurements were made 
with the transmitter and receiver placed a few feet apart with the local 
laser in the receiver turned off, so that only noise contributions from 
the transmitter laser, detector and receiver electronics would be 
present. The system was then operated in the same manner over a 100 
meter path through an enclosed tunnel one meter in diameter. The noise 




Figure 7. Block Diagram of Data Acquisition System. 
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of the local oscillator laser was determined by operating it into the 
detector in the absence of an incoming signal. The noise in all cases 
was found to be sufficiently low so that it would be negligible compared 
to the expected variation due to atmospheric scintillation. The linearity 
of the mercury cadmium telluride detector was determined by noting 
changes in DC voltage output for different power levels. The laser 
power was measured on one side of the, germanium beam splitter using a 
Coherent Radiation Laboratories power meter and the DC variation in the 
detector was- measured by a digital volfmeter. for incident power levels 
less than 300 mw. , the detector output was found to be linear. 

Alignment of the system over the 3.2 km. path proved to be a 
difficult task. Our first attempt was to bore sight a 60 power tele- 
scope mounted on the transmitter case, with the output beam. The idea 
was to aim the transmitting unit on the mountain at the laboratory 
periscope. This method worked over the 100 meter tunnel quite well, but 
the bore sight became misaligned when the transmitter was transported to 
the mountain. The second method employed to align the system involved 
the use of two visible lasers. The transmitter unit, now mounted on the 
observatory telescope stand on Madkin Mountain, was aligned by placing 
an argon laser directly in place of the receiving unit in the laboratory. 
The bright beam could easily be detected by the eye at the transmitting 
terminal. The position of the argon laser was adjusted until its beam 
was intercepted by the objective of the transmitter unit. The optical 
system of the transmitter was then adjusted so that the visible laser 
was focused onto the output aperture at the transmitter laser. Using 
heat sensitive paper as a position indicator for the infrared beam, the 
visible light and the invisible beam were made to coincide at two points 
in the optical system. Alignment of the receiving unit was accomplished 
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in a similar manner except that a visible laser could not be mounted in 
the same position as the now aligned transmitter unit. A small helium 
neon laser was mounted with a telescope on a tripod and located as near 
to the transmitter as possible. The telescope was then bore' sighted to 
the helium neon laser. The laser-telescope arrangement was pointed by 
locating the top periscope mirror at the laboratory. A corner re- 
flector located at this mirror enabled a more precise aim as the re- 
flection of the red light could be seen with the telescope. The 
visible beam was then focused onto the detector in the receiver. With 
minor adjustments to the transmitter laser mount, close alignment was 
attained for the system. 

C. Measurement Procedure 

Scintillation measurements were made with the receiver laser 
inoperative so that only light from the transmitter and from the sun's 
reflection off the observatory dome were intercepted by the receiver. 

This background light was cause for great concern since accurate 
measurements of the variations of the laser light could not be made in 
its presence. Since we could not filter out this unwanted light, it was 
necessary to record it. To accomplish this scintillation measurements 
were made by chopping the transmitter beam at 90 Hz by means of the 
mechanical chopper located in the transmitter. Figure 8 shows a sample 
of this chopped signal. 

Scintillation measurements were made in 90 second segments. A written 
log was kept on each data run giving the date and time of day, analog tape 
number, the channel number, aperture size, and weather conditions such ' 
as temperature, wind velocity and humidity. The intensity signal for 
each data run was read on the analog tape into marked time slots. One 
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of the channels was marked every 90 seconds with an audio tone code so 

1 

that approximately 33 slots were available on each data channel. When 
sufficient signal was being received, data runs were made in succession 
using five different size apertures. The aperture size varied from two 
to ten cm. in diameter. The runs were made in succession to insure that 
atmospheric conditions remained constant over each set of five runs. 

Since the signal level for small aperture sizes was below that for 
the larger apertures, it x<ras necessary to adjust the AC amplifier from 
one run to another using the monitor oscilloscope. Also, it was 
important to keep the dewar on the detector filled with liquid nitrogen 
so that the detector temperature would remain constant. 

, Measurements of the signal to noise ratio and distortion were made 
on the communication system operating over the optical range. The si gn al - 
to-noise ratio of the heterodyne action was measured at the receiver by 
extracting the 10 Mhz beat note between the received signal and the local 
oscillator after it had passed through one stage of amplification. The 
signal was detected with a simple diode circuit and the resulting 
voltage recorded by the data acquisition system. These measurements were 
also made for different aperture sizes. 



CHAPTER IV 


DATA REDUCTION 

A. Analog-to-Digital Conversion 

All' measurements were recorded on analog tape. In order to process 
this data it was necessary to convert it to a digital form suitable for 
computer input. The very large quantity of data collected necessitated 
electronic conversion to digital magnetic tape which could then be read 
. directly into the computer. 

The first step in the conversion from analog- to-digital form con- 
sisted of recording a timing signal on a reserved channel of the analog 
tape. Eight standard time signals are broadcast by Marshall Space Flight 
Center's Computation Laboratory on a frequency of 226.5 Mhz. The time 
signals are subcarrier multiplexed with pilot frequencies between 2.3 KHz 
and 70.0 KHz. The signal chosen was a rectangular pulse with a one- 
second repetition rate broadcast on a subcarrier frequency of 52.5 KHz. 
This signal is designated as 100/1000 AMR-D-5, and is coded with the 
time of day in Greenwich Mean Time by pulse width modulation. This timing 
signal had no relation to the actual time the data was recorded. 

The analog tapes were then read into a cathode-ray-oscillograph. 

Four data channels, the timing channel and the marker channel were re- 
corded simultaneously. The oscillograms were inspected and the sections 
of the tape to be digitized were selected. At this time bad data was - 
identified and eliminated. The starting and stopping time for each time 
, interval selected was read from the timing channel and recorded on the 
Computation Laboratories instruction forms . The analog-to-digital con- 
verter system was set to convert only those time intervals selected by 
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reading the timing channel. Approximately 30 one-minute segments were 
selected from each channel. As ■ the - ' running time for the analog tape 
was about 45 minutes , about two-thirds of the total data recorded on a 
tape channel was converted. 

The analog-to-digital conversion was performed by Marshall Space 
Flight Center's Computation Laboratory using an Astrodata type con- 
verter. Five data channels were digitized simultaneously. The five 
channels were sampled alternately at a rate of 5000 samples per second' 
which amounted to sampling each channel at 1000 samples per second. The 
resulting binary coding was recorded on seven channel digital tape in a 
multiplexed format . 

The digital tape format consisted of a ten bit binary word so that 
2048 levels were available to represent analog signal levels between —.5 
volts and +.5 volts. Since seven track tape was used, each sample re- 
quired two tape characters consisting of a ten bit word plus the sign 
bit, blank bit and two parity bits. The data was recorded in records of 
2004 tape characters. The first four tape characters contained the time 
as read from the timing channel at which the first sample was taken. The 
remaining 2000 characters contained 1000 sample points, 200 from each 
channel, alternating between channels. Subsequent records were written 
until the segment was completed, at which time an end of file mark was 
written; therefore, each- file on the digital tape contained one time 
slice or data run from the analog tape. In addition, the first record 
of each file was an identification record of 24 tape characters which con- 
tained the tape number and other information. 

B. Description of Computer Program 

A program was written for an IBM- 360-50 computer to reduce the data 
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stored on magnetic tapes- This program reads data from the tape, changes 
the binary format to a fortran compatible form, then calls its various 
subroutines to perform the analysis. The principle problems which were 
encountered in writing this program concerned formatting the data for 
the computer and extracting the actual light intensity signal from the 
modulated square wave. When the signal was extracted it was either 
stored for spectral analysis or is used to construct a histogram. A 
Fourier transform subroutine or statistical subroutine was then called 
to analyze the signal. 

Development of a routine to extract the desired signal proved to 

be somewhat difficult since the sampling rate during digitization could 

not be accurately synchronized with the period of the square wave. The 

sampling rate of 1 KHz and the chopping rate of 90 Hz should yield 
* 

approximately 10 samples per cycle of the square wave. In actuality, 
the number of samples per cycle varied between ten and twelve due to the 
sampling rate not being an integral multiple of the square wave frequency. 
The routine was designed to determine whether a particular data sample 
was a base point (i.e., from the part of the square wave when the laser 
beam was blocked by the chopper) or a signal point (when power was being 
received from the laser beam) . The problem was further complicated by 
the fact that the rise and fall times of the square wave were non- 
negligible so that about one percent of the data points were sampled 
during the switching transient and should be neglected. In addition, 
some of the data contained an occasional noise spike which should be 
eliminated. It was decided that the elimination of these spikes would 
not adversely effect the validity of the analysis so provisions for 
eliminating them were also included in the program. 
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The extraction routine operates basically as follows: A pre- 

processing routine reads the data from the magnetic tape and stores a 
record containing 200 data points into a common array. To begin the 
analysis twenty data points from the array are selected and their maximum 
and minimum value computed. Two limits, and are then set by the 
relations 
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P, (A 
1 max 
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where A and A . are the maximum and minimum values of the first 
max mm 

twenty points, and P^ and are constants between zero and one half. 
Since the signal was inverted when it was recorded on the analog tape, 
the base line is greater than the signal, hence a particular point 
greater than is considered a base point, if it is less than L 2 it is 
considered a signal point. Points lying between and are assumed 
to be from the transient portion of the wave form and are neglected. The 
routine takes each point successively and determines if it is a base 
point, a signal point or neither. As a preliminary to processing, the 
first twenty points are scanned and the beginning of a base line segment 
of the wave form is found. Then new limits are set on the next 15 points 
and they are scanned and grouped into three arrays, a base line segment, 
a signal segment and a second base line segment. Each array may contain 
up to ten points . At this time another routine is used to compute the 
signal amplitude of the square wave for the group of signal points (as 
will be described later) . The second group of base points is transferred 
into the first array, new limits are set using the next ten data points. 
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and a new group jof signal points- and base points are found to fill the 
second and third arrays, and finally their amplitudes computed. This 
process is continued until the 200 points from the first record have been 
used. At this time the routine pauses while the next record is read in. 
Processing then continues until, the number of records called for (up to 
300) have been processed. 

This routine also contains several checks to handle possible 
irregularities in the data. During the search for either base or signal 
points if more than ten consecutive points are found, the routine will 
request that the next record be read in, which means the remainder of 
the bad record is discarded. Also if the number of unsuccessful scans 
while in the base or signal search phase exceeds ten, the routine will 
enter an error recycle phase. In this phase the routine skips 20 points 
and resumes processing as if it were at the beginning of a new run. If 
the error recycle phase is entered five times in a given record, a re- 
quest for the next record will be executed. When a new record is re- 
quested due to an irregularity in- the record being processed the routine 
again treats it as it does the first record’ of a new run. In addition 
to the above, if three or less base or signal points are found in a 
given search, the error recycle phase is entered. 

During processing, a record is kept of each irregularity encountered 
and this information is printed in tabular form when the processing of 
a run is completed. If an excessive number of irregularities occurs in 
a given run, the results of that run must be suspected. 

The three arrays containing base and signal points found by the 
extraction routine are passed into a routine that determines the actual 
amplitude of the signal. Three methods for computing the amplitude were 
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tried: The first method took to base line for a- group of signal points 

as the average of all the points in the group of base points preceding 
it and the ones following it. This is equivalent to considering the 
background light during the time when the laser beam intensity was being 
recorded to be the average of the -background light recorded during the 
half-cycle immediately preceding the signal and the half-cycle immediately 
following it. The difference between the signal point and the average 
base line was taken as the amplitude of the laser beam at that instance. 

The second method of amplitude calculation considered was to re- 
construct the base line by fitting a least-mean-square curve to the base 
points . This method produced erratic results and -also required additional 
computer time and was abandoned. 

In the third approach, the difference in the first signal point in 
a group and the last base point preceding it is taken as the signal 
amplitude. The difference in the last signal point in the group and the 
first base point following it gives a second amplitude. This method 
yields only two amplitudes per cycle but has the advantage that they are 
evenly spaced. Since the difference in time between signal and base 
points is very small, this method eliminates the problem of the unknown 
base line over the signal interval. 

The computer program contained both the first and third methods of 
signal amplitudes calculations. The method to be used was selected by 
a "parameter read during execution. The values computed by this routine 
were either stored in an array to be used in the spectral analysis, or 
used to construct a histogram. Histograms are sometimes referred to as 
being the probability density function for discrete variables. 

The final segment of the program consists of the analysis routines. 
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The statistical routines accept the histogram for the intensity 
fluctuation which has been generated and computes the scintillation 
statistics. The routine computes and lists the corresponding value of 
the log-amplitude as defined by 

L = \ in [I /T] 4-3 

where L and are the log-amplitude and the intensity for the ith class 
interval and I is the mean intensity. The frequency for each class and 
the cumulative probability are also listed. The routine also calculates 
the mean, standard deviation, skewness, and kurtosis for both the in- 
tensity and log-amplitude distribution. A chi-square test that checks 
the intensity distribution for a normal fit and the log-amplitude for a 
log-normal fit is also included. Appendix A includes a typical computer 
print out of the .statistical analysis. 

The program includes routines to perform a spectral analysis on the 
intensity scintillation. In calculating the scintillation spectrum 

g 

2 (8192) values of the beam intensity were extracted from the raw data 

using the preprocessing routines described previously. These values 

represent the intensity fluctuations of the beam sampled at 180 Hz rate. 

The spectrum of the sampled data was computed using the Cooley-Tukey 
9 

algorithm for fast Fourier transforms. The resulting spectra were dis- 
played by means of a plot routine. In cases where data points are dis- 
carded due to irregularities in the extraction routine, the omitted points 
were replaced with zero values. The purpose of this was to preserve the 
phase relationship of the signal, which is important in the integration 
operation of the Fourier transform. 
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C. Program Check and Parameter Adjustment 

To test the program several records were read from a magnetic tape 
and printed out for inspection. Each sample was classified either as a 
signal point, or a base 'point, or as a point from the transient portion 
of the waveform. This classification was purely subjective, yet in 
inspecting the data there was usually no questions as to how a point 
should be categorized. The Same data was then read into the computer. 

A print out of all the pertinent variables at each decision branch of the 
extraction routine was obtained. The program was then run several times 
varying the parameters P^ and P^ in equations 4-1 and 4-2 between runs 
and the results compared with the subjective analysis. Data having 
irregularities was also processed and the results compared to determine 
whether or not a segment should be omitted. Using these comparisons as. 
a criteria,- the parameters P^ and P^ were set at 0.05 and 0.10 respectively, 
Therefore a point within 5% of the maximum base point or 10% of the 
minimum signal point would be retained while points between these limits 
were discarded. 

The statistical routines were checked out with sets of numbers which 
had a known probability distribution function. A routine readily available 
in the IBM Scientific Subroutine Package for the IBM 360 computer was used 
to generate a large set of normally distributed numbers . These numbers 
were then processed by the statistics routines in the same manner as the 
intensity fluctions were to be processed. The part of "the routine which 
tests for the normal distribution fit gave very positive results, and 
that for the log-normal fit gave negative results. The test routine was, 
then altered to generate numbers with a log-normal distribution. The 
results of the statistical analysis of this data strongly indicated a 
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log-normal fit. A description of the log-normal generator is given 
in Appendix B. 

D. Atmospheric Structure Constant and Aperture. Averaging 

An important parameter in the statistical model for atmospheric 

studies used in current literature is the refractive index structure 
2 

constant C , defined as a constant of proportionality in the relation 

D n ■ <^tn 1 (r 1 )-n 2 <r 2 )]^> = C 2 r 2/3 

r = r 2 - r 1 4-4 

where D is called the structure function and is a measure of the 
n 

deviation of the index of refraction at two points. separated by a 
distance r. is actually a measure of the strength of the turbulence 

Fried ^ gives a relation involving for an infinite spherical wave 

propagating a distance z in a turbulent atmosphere. The relation is 

C.(0) = .124 K 7/6 Z 11 ^ 6 C 2 ’ 4-5 

Jt n 

where 0^(0) is the log-amplitude variance. Equation 4-5 can be used 

2 . 

directly to obtain by noting that the standard deviation of the 
log-amplitude distribution which we have computed is the square root of 

V 0) - 

A finite receiving aperture has the effect of averaging the inr- 
tensity fluctuating from various parts of the wave front thereby re- 
ducing the variance of the scintillation. 

The effect of aperture averaging can be allowed for by using re- 

11 12 

lations developed by Fried ’ viz. 
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^ s 2 = rf D 2 ] 2 0 c i(°) 4-6 

Where a g is the signal variance which corresponds to the square of the 
standard deviation of the intensity fluctuation, D is the diameter of 
the receiving aperture, and 0 is an aperture averaging factor given by: 


0 
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H(p/D) is the optical transfer function of a circular aperture 

H(p/D) = cos 1 (p/D) - p/D [l-Cp/D) 2 ] 1 ^ 2 

and C^(p) is log-amplitude covariance given by: 
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In the last expression a^ and are the expansion coefficients 
for the modified confluent hypergeometric function and are given by 
the recursion relations 


a n = -a n _ 1 { (2n - 23/6) (2n - 17/6) / (2n-l) (2n) } 

b n - -b n _ 1 { (2n - 11/6) (2n - 17/6) (2n-l)/(2n) (2n+l) 2 } 
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where 

a = 1 b = 6.84209 

o o 

The intensity variance C^(0) can be related to the log-variance by: 

C I (0) = I q 2 [exp(4C £ (0)> - 1] 4-11 

' 2 2 

Equation 4-11 specifies C„(0) in terms of a and I . I and <r are 

l s o - o s 

the mean and standard deviation of the recorded intensity. Combining the 
above equations we have: 
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where f(kp/4z) is the summation given in equation 4-9. Since o^/I is 

an experimentally determined constant, equation 4-12 is an integral 

equation for (0) . A special computer program was written to solve 

equation 4-12. The technique used is to evaluate the integral in 

equation 4 J -12 for a number of trial values of C (0) using a fourth order 

a 

Runga-Kutta integration. This gives a table of —■ as a function of ' 

o 

C (0) . From this table the value of C. (0) corresponding to the measured 
* a £ 

g 

— is determined using Lagrange-Hermite interpolation formula, 
o 


value of 



CHAPTER V 


RESULTS 

A. Probability Density Eunction for Intensity Scintillation 

It has been customary in the literature to test the hypothesis of 
log-normality of scintillation data by plotting the cumulative prob- 
ability function of the log-amplitude against a "probability scale" such 
that if the data is log-normal the resulting curve will be a straight 
line. The same method can be applied when testing the intensity amplitude 
for a normal distribution. Since this test would require considerable 
time if it were applied to every run, it was necessary to use a test 
that could be incorporated into a computer program in an efficient manner 
in order to determine which distribution each run more closely fit. 

The necessary statistical parameters to make this test were calcu- 
lated as described in Chapter III and were available on punched data 
cards. The skewness coefficient was chosen as the parameter to indicate 
the type of distribution. The skewness coefficient will ideally have 
zero value for perfectly normal or log-normally distributed data; however, 
for real data we expected a small value but somewhat greater than zero. 

We have chosen the skewness coefficient as the measure of closeness 
of fit in preference to the chi-square test since the chi-square depends 
upon the number of class intervals in the. sample while the skewness 
does not. To use the chi-square on a comparative basis would require the 
generation of a table of chi-square for all possible number of class in- 
tervals and would lend to undesirable complexities in the computer pro- 
gram. 
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In order to use the skewness as a criteria for categorizing the 
data it was necessary to set bounds on its value. Since the lower 
bound is zero, it was only necessary to determine the largest value that 
the skewness could attain which would represent a suitable fit. This 
was accomplished by plotting several graphs of the cumulative probability 
for runs with values of skewness ranging from .02 to .79. The chi- 
square test made on each run was examined to insure that the overall 
results of the analysis were consistent. A selected number of these graphs 
are given in Figures 9-17. An inspection of these plots will show that 
the cumulative probability curve becomes nonlinear with increasing 
values of skewness . An inspection of many such plots indicated that the 
curve deviated from linearity much faster when values of skewness be- 
came greater than 0.15. For values from 0.0 to 0.15 the curve remained 
approximately linear. 

* 

The results of the statistical analysis of each run was categorized 
in the following manner: If the skewness for the log-amplitude data is 

less than or equal to 0 ..15 , the run was categorized as being log- 
normally distributed. If the value of the skewness meets the above re- 
quirement for the intensity data the run was considered to be normally 
distributed. If both skewness coefficients were greater than 0.15, the 
run was put in a "neither" category. In addition we have had to in- 
clude a category for those distributions which were sufficiently close 
to both a normal and a log-normal distribution that we could not 
distinguish between them. These runs which have approximately the same 
skewness for both distributions have been designated as "both". 

The results of the categorization for 196 runs are shown in Table 
I. Only runs with a small number of preprocessing irregularities have 
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AMPLITUDE 
(arbitrary scale) 


99.8 99.9 99.99 



APERTURE 

SIZE 

(CM) 

NUMBER OF RUNS 

TOTAL 

LOG 

NORMAL 

NORMAL 

NEITHER 

BOTH 

2 

4 

2 

2 

0 

0 

4 

17 

13 

; 2 

1 

1 

6 

30 

10 

6 

(3 

1 

8 

43 

1 

1 2 

14 

. 1! 

6 

10 

102 

3 3 

2 3 

37 
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Table 1 , Results of Statistical Analysis of Scintillation Measurements. 
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been tabulated. From this table we see that the two cm. aperture runs 
were divided equally between normal and log-normal distributions. The 
results for this aperture size are not considered significant since there 
were only four .runs with small irregularities, and these were taken with 
a very weak signal. The four cm. runs, which had a reasonably small 
aperture, yet the recorded signal was strong, show a strong tendency 
toward log-normality. As the aperture increases the number of runs which 
differ from log-normal also increases. This behavior may be due to the 
effect of aperture averaging. Since this process is additive in nature, 
it would cause the distribution to tend toward a normal curve according 
to the Central Limit Theorem. This combined with the difficulty in 
distinguishing normal data from log-normal data could lead to results of 
the type exhibited by the larger aperture runs. 

B. Calculation of Atmospheric Structure Constant 

The refractive index structure constant was calculated using 
equation 4-5 and taking the measured value of the log-amplitude variance 
C (0). This equation does not allow for aperture averaging effects. 
Values of C^ obtained ranged between 1.6 x 10~ 8 M 3 ^ 3 and 8.7 x 10~ 8 M 3 ^ 3 
for a group of runs where the aperture was varied rapidly from two to 
ten cm. The structure constant computed from two cm. aperture data was 
usually larger than that computed from ten cm. aperture data by a factor 
of from two to four. The time required to collect the data for such a 
group of runs was less than ten minutes. As a comparison, several 
groups of runs were made in a similar time period holding the aperture 
constant. The observed variation in the structure constant for these 
rims was usually about ten percent. This clearly indicates that the 
observed variation in the structure constant was due to aperture 
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averaging effects and not to changes in the strength of turbulence be- 
tween runs. 

Using the techniques described in Chapter III, section D, calcula- 
tions for the structure constant have been refined, allowing for the 
effects of aperture averaging. For the 212 segments of data which were 

analyzed, the values of the corrected structure constant ranged from 
-7 1/3 -7 1/3 

5.8 x 10 M to 9.0 x 10 M . These values are characteristic of 
very strong atmospheric turbulence, which agree with our subjective 
observations of the scintillation while the data was being taken. 

Table II shows four typical sets of runs for various aperture 
sizes. As can be seen, the variation in the structure constant is 
significantly reduced when the effects of aperture averaging are in- 
cluded. 

It should be noted that the inclusion of aperture averaging effects 

has a tendency to overcorrect the structure constant variation. This 

could be an indication of a systematic error in recording the data. 

Another possibility is that the deep scintillation conditions under which 

the experimental data was collected produced saturation effects which have 

not been considered. On the other hand the validity of the basic approach 

to the problem of aperture averaging in terms of the structure function 

13 

has been questioned . Therefore it is possible that the aperture 
correction we have used is not valid. In any case, this method of com- 
pensating for aperture averaging effects is a good approximation since 
the structure constant is far more consistent when corrected than when 
aperture averaging effects are neglected. 

C. Scintillation Frequency Spectrum 

The frequency spectrum has been computed using the Fast Fourier 
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Transform techniques described above. This calculation was performed 
only on selected data runs due to the time required for such a calcula- 
tion. Also the computed spectra were very consistent so it was felt 
that analysis of additional runs would yield little additional informa- 
tion. • 

The computed spectra cover a range of frequencies from DC to 90 Hz, 
the upper limit being set by the 180 sample per second sampling rate. 

A typical spectrum is shown in Figure 18. Although the computer plotted 
frequency components out to 90 Hz only components out to 20 Hz are shown 
in the figure for the sake of clarity. Above 20 Hz the spectrum con- 
tinued to decrease linearly so that no appreciable frequency components 
above 40 Hz were observed. 

D. Heterodyne Detection 

The effect of atmospheric turbulence on the performance of the 
equipment operating as a heterodyne communication system was. investigated. 
This was accomplished by recording the amplitude of the 10 Mhz 
heterodyne signal at the output of the I. F. amplifier. Also the trans- 
mitter was modulated with a 1 Khz signal which was. recorded at the out- 
put of the receiver's F. M. descriminator . It was found that neither 
signal showed any effect attributable to atmospheric turbulence large 
enough to be accurately measured. Even under conditions of deep 
scintillation encountered during the course of this experiment, the 
atmospherically induced noise was of the same magnitude or smaller than 
system noise. It was also found that clear voice communications were 
possible over this range under the worst conditions of scintillation 
encountered. 

While these results clearly indicate the feasibility of using a 
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heterodyne CO^ laser system to communicate through the earth atmosphere, 
they were somewhat disappointing since they did not permit a quantitative 
measurement of noise induced by atmospheric turbulence. 





CHAPTER VI 


SUMMARY AND CONCLUSION 

The results of the scintillation measurements made on the 10.6 
micron wave.length laser beam tend to confirm the log-normal model for 
small receiver apertures. The data for the larger apertures did not 
seem to fit the log-normal or normal models with any consistency. One 
possible explanation for this could be aperture averaging. It is 
possible that the larger aperture sizes were not large enough with re- 
spect to the correlation distance to cause the distribution to be normal,' 
but yet large enough to cause the distribution to deviate from log- 
normally. This in addition to the difficulty in distinguishing between 
the two distributions could have caused the results to be inconclusive. 

The value of the refractive index structure constant computed was 
found to lie within the range of values for this constant as calculated 
by Fried. The value of this constant was found to decrease with in- 
creasing aperture size. Equations developed by Fried were used to correct 
the structure constant for aperture effects. This technique seemed to 
give a slightly larger value of the aperture averaging effect than was 
observed. Although this may indicate an inaccuracy in the theoretical 
expression, a systematic error in the experimental measurements cannot 
be firmly ruled out. These calculations were significant in that they 
indicated in a quantitative manner the nature of aperture averaging. 

The spectral analysis indicated that low frequency components of 
the scintillation were predominant. The magnitude of the scintillation 
decreases linearly with increasing frequency. Above 20 Hz the 
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scintillation is negligible. 

The feasibility of optical heterodyne communication at 10 microns 
through extreme turbulence was demonstrated. To our knowledge, this was 
the first system of this kind to be operated through the atmosphere at 
this path length. We feel that the successful operation of the 
communications system under extreme scintillation conditions was a 
significant result. 

It would be a great advantage in this type experiment to reduce the 
data immediately after it is taken. Information gained from the speedy 
reduction should give the experimenter a knowledge of how the system is 
performing and aid in making better measurements . 

As a result of performing this experiment and surveying the results 
of other investigations, it is clear that the atmospheric problem is far 
from being completely solved. The log-normal model needs to be further 
verified for other wave lengths. The variance and the structure constant 
should be investigated under as wide a variety of weather conditions as 
possible and should be correlated to the variations of the meterological 
parameters . 

After having carried out this type experiment, the need for several 
refinements in the procedure was realized. Before any data is taken, a 
thorough analysis of system noise should be made. Sensitive noise 
measuring instruments should be employed. The noise of the transmitting 
laser should be recorded simultaneously with scintillation measurements. 
The background light effect should be further studied and perhaps a 
method other than signal chopping used to handle the problem. Mechanical 
stability of equipment is a problem that needs investigation since even 
very small vibrations could cause the beam to shift out of alignment 
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with the receiver. 

Aperture averaging .effects need to be investigated with many 
different size apertures for all wave lengths. The relationship between 
correlation distance and aperture size needs to be determined. The 
determination of correlation distance in itself would be an interesting 
experiment. 
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APPENDIX A 


Introduction. 

This appendix further describes the computer program employed to' 
reduce the atmospheric data. The program is written in Fortran IV 
language for the IBM— 360 Model 50 computer at the University of Alabama. 

The program operates in the following manner: The MAIN or supervisory 

routine accepts instructional and operating data for the program. It 
reads the atmospheric data from the magnetic tape and calls subroutine 

EXT to extract the intensity signal. EXT' calls on subroutines LIMIT, 

< 

FILL, HIST, and AMPX to perform the extration. When a file of data has 
been processed MAIN calls subroutine PRINT to print a table of the 
irregularities that occurred during extraction. MAIN then calls sub- 
routine STAT and/or FFT to perform the statistical and spectral analysis. 

STAT calls on subroutine CHI to perform the chi-square test which in turn 

) 

calls on a Simpsons rule integration subroutine SIMP. FFT calls on two 
package subroutines FOURT and PLOT which perform and plot the spectral 
analysis . 

Routine Main 

Routine MAIN's first step is to read the identification record from 
the magnetic tape. This is done by calling subroutine RID. RID 
actually does the reading and stores the data in a common array to be 
printed by MAIN in the next step. The- numerical instructions and operating 
constants are then read in as data on punched cards . These are read in 
as variable names and are used in the various subroutines and are dis- 
cussed as part of the description of these subroutines. 
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Processing actually begins when MAIN reads in 4 additional in- 
structional constants on punched cards. These are given according 'to 
their variable name as 

NAME FUNCTION 

FILES — The number of the data file to be processed. 

ISTART - The record number within the file from which data 
will start being taken. 

ICHNO - The channel number of the magnetic tape to be 
processed. 

IRCEND - The record number at which processing is to stop. 

The values of these numerical constants enable the user to process any 
record or segment of records within a data file on any of the five 
channels . Each time a new file is to be processed these variables must 
be read in for the new file. The read statement for these variables is 
in a loop so that when processing of a file is completed the program re- 
turns to this statement to get instructions for processing the next file. 
After the last file has been processed the program is stopped by 
entering a negative number for the variable FILES. 

The program now moves the tape to the desired file and record by 
calling subroutine REDREC. When the first record to be processed is 
found, it is stored in a common array IBUF. MAIN then transfers the 
data in IBUF into a work array P. REDREC is called again to read the 
next record which is transferred by MAIN into an auxilary array AUX. 

This is done in order to have the next record available in an array 
since it is sometimes necessary for the program to' "look" into the 
following record before processing in' a record is completed. 

The record of data in array P is now processed by calling sub- 
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FORTRAN IV G LEVEL 1, MOD 4 MAIN 


DATE = 70122 02/35/57 


0001 

0002 


0003 

0004 

0005 

0006 

0007 

0008 
0009 
00 10 
0011 
0012 
0013 
■ 0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 
0023 


0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 


INTEGER FILES, FUCK 

COMMON P (1000), AtJX (1000) ,INDO(10) ,NCNT(300) ,JOVF,JUNF, 

♦BASE (2, 10) , SIG (10) ,IBLOC{2,10) , KK, HM,XL90 ,XL1 0 
*, DUMMY (2700) 

COM MON/A A A/IPRINT 
COMMON/AP/AHP (10000) , NDATA 
COMMON/RLC/IRTN 
COMMON/RAT/IRCEND 
■ COM MON/GO/N 1,N3,N4»FM,IDO 
COMMON/NPTS8/N8 
COM HON/PB/IE (7,300) 

COMHON/UUU/NPNCH 
DIMENSION IBUF(IOOO) 

DIMENSION IA (3) ,IB (5) 

CALL RID (11,12, IB) 

IF (II) 155,150,155 
155 PRINT 1 57 ,FILCK 

157 FORMAT ( 1 *** END OF FILE ♦***,18) 

150 IF-CI2) 158,161,158 

158 PRINT 159 

159 FORMAT {' *** READ ERROR *♦**) 

161 CONTINUE 

NTAPE=IB (1) 

PRINT 160, IB 

160 FORMAT (• TAPE NUMBER*' ,1A4, 

1* IDCAL= * , 1A4, 

1* BASE 55 * , 1A4', 

1* NOMBER OF CBANNELS=* , 1 A4, 

1* NUMBER OF SCANS*', 1A4) 

FILCK= 1 

READ (5, 9921) FHAX,DT,FM, N1 ,N3, N4,IDO,ISSB 
9921 FORMAT(3B10. 0,5110) < 

READ (5, 3350) IPRINT, NPNCH 
3350 FORMAT (2110) 

C READ INSTRUCTIONS AND TOLERANCES • 

READ (5,36 0) IBTYP,IFLAG, NCI, NOSCAN, NOCHAN, LI, L2 

360 FORMAT(7I 1 0) 

READ (5, 1991) TOL1 , T0L2, TOL3 
.1991 FORMAT (3E10. 0 ) 

PRINT 361 ,IBTYP,IFLAG, NCI, NOSCAN, NOCHAN, TOL1 , TOL2,TOL3 , Li, L2 

361 FORMAT ( * TYPE BASE CALCULATION- ********************************* 

****** ***** IBTYP ***** *110/, 

*■ TYPE STATISTICS CALLED FOR ********************************* 
****** IFLAG ***** ',110/, 

*' NUMBER OF CLASS INTERVALS ********************************** 
****** HCI ***** » ,110/, 

*» NUMBER OF SCANS PER RECORD ********************************* 
****** NOSCAN ***** ’,110/,. 

*' NUMBER OF CHANNELS ON TAPE ********************************* 
****** NOCHAN ***** ',110/, 

** TOLERANCE ON BASE LIMITS *********************************** 
****** TOL 1 ***** ’ , FI 0. 3/, 

*' TOLERANCE ON SIGNAL LIMITS ********************************* 
****** TOL2 ***** ' ,F10.3/, 

*> NUMBER OF SIGNAL AND BASE POINTS REQUIRED PER CYCLE ******** 
****** TOL3 ***** * ,F10-3/, 
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FORT RAH 


0035 

0036 


0037 

0038 

0039 
004,0 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 


0053 

0054 

0055 

0056 

0057 

0058 

0059 

0060 
0061 

0062 

0063 

0064 


IV G LEVEL 1, MOD 4 


MAIN 


DATE = 

70122 

02/35/57 

*> NUMBER 

OF 

PASSES REQUIRED TO ABORT 

CYCLE 

******************* 

****** xi 


***** 

*,no/. 




*’ NUMBER 

OF 

EASE POINTS 

REQUIRED FOR 

SIGNAL 

POINT SEARCH 

***** 

****** x2 


***** 

>,110) 





PRINT 911 1,NPNCH,IPRINT,ISSB f FMAX,DT,FH,N1 ,N3,N4,IDO 
9111 FOHMATf PUNCHED OUTPUT FLAG **********************.************* 
****** ***** NPNCH ***** ’,110/, 

*» PRINT ERROR TABLE FLAG ************************************* 
****** iprint ***** ',110/, 

* • TYPE ANALYSIS ROUTINE RANTED ******************************* 
****** X5SB ***** *,110/, 

*t INFORM ATIOH FOR FFT ROUTINE ******************************** 
****** FM ax ***** ' ,F1 0-3/, 65X, 1 ***** DT ***** 

** ,F10. 3/, 65X, ****** Ffl ***** * , F 10 . 3/, 65X, 1 ***** N1 

* ***** 1 ,110/, 65X, ****** N3 ***** 1 ,1 10/,65X, ' **. 

**** N4 ***** 1 ,110/, 65X, • ***** IDO ***** *,I1 

* 0 /) 

C READ DATA FOR TAPE PROCESSING 

LFILE=0 

51 1 READ (5,882) FILES, ISTART,ICHNO. IRCEND 
N8=0 

NDATA=0 . 0 
DO 8892 EJ=1 , 7 
DO 8892 JI=1, 300 
8892 IE(IJ,JI)=0 

IRECNI=ISTART-1 
IEX IT=0 

IF (FILES- LFILE) 22,22,987 
987 IREC=0 
22 CONTINUE 

LFILE=FILES 
882 FORMAT (413) 

PRINT 901 , FILES, TCHNO,ISTART, IRCEND 
901 FOR BAT ( 1 H 1 , 1 PROCESSING FILE **********>,15/,, 

*< CHANNEL *****************************< ,15/, 

*• BEGIN AT RECORD *********************> ,15/, 

*' STOP PROCESSING AT RECORD ***********>, 15) 

IRTN=1 

IF ( FILES) 991,9,9 

991 PRINT 992 

992 FORMAT ( 1 0X, * ****** PROGRAM STOP *♦***>) 

STOP 

9 CALL REDREC (FILES, IREC,FILCK,IBUF, TIME, N, NOSCAN, NOCHAN, ICHNO) 

C MOVE TAPE TO DESIRED FILE 

IF (FILCK-FILES) 9,887,887 
C READ RECORD, STORE RECORD IN ARRAY IBUF 

887 CALL REDREC (FILES, IREC,FILCK, IBUF, TI ME, N, NOSCAN, NOCHAN, ICHNO) 

C HOVE TAPE TO DESIRED RECORD 

IF ( IREC— ISTBRT) 887,889,889 
C STORE CONTENTS OF IBUF INTO MAIN ARRAY P 
889 DO 888 L=1,N 

888 P(L)=IBUF(L) 

C READ NEXT RECORD INTO IBUF 

CALL RED REC (FILES, IREC, FILCK, IBUF, TIME, N, NOSCAN, NOCHAH, ICHNO) 

C STORE CONTENTS OF IBUF INTO AUXILIARY ARRAY AUX 

DO 500 L=1,N 


0065 
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0066 

500 

AOX (L) =IBUF (L) 

0067 


GO TO 501 

0068 

502 

DO 503 L=1,N 

0069 

503 

P {L)=AOX (L) 

0070 


CALL RED REC(FILES, rREC,FILCK,IBUF,TIHE,N, NOSCAN, NOCHAN, ICHNO) 


C 

STORE IBOF 1 8 TO AIJX 

0071 


DO 504 L=1,N 

0072 

504 

AUX (L)=IBUF (L) 

0073 

501 

CALL EXT (N ,L 1 ,L2,IRECNT, TOL1 ,TOL2,TGL3 ,NCI,IBTYP) 


C 

STOP PROCESSING ON DESIBED RECORD 

0074 


IE(IEXIT) 1000,1000,507 

0075 

1000 

IF {NDATA.EQ. 10000) GO TO 507 

0076 

661 

IF(IREC-IRCEND) 502,221,221 

0077 

221 

DO 2293 1=1, N 

0078 

2293 

P (I) =ABX (I) 

0079 


IEXII= 1 

0080 


GO TO 501 


C 

PRINT ERROR TABLE 

0081 

507 

CALL PRINT (1) 

0082 


IF (ISSB.EQ.2) GO TO 2514 

0083 


CALL FFT (N , DT ,FHAX) 

0084 


IF(ISSB.EQ.I) GO TO 511 

0085 

2514 

CALL ST AT.(NTAPE. ICHNO. FILES. NCI. IFLAG) 

0086 


GO TO 511 

0087 


END 
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routine EXT. When EXT completes the processing it returns program 
control to MAIN. The program checks the variable IEXIT to see if the 
record just processed is the last record in the file. If it is not, 
the program checks to see if the next record is the last. If this re- 
cord is not the last, the program calls in the next record and continues 
processing. However, if it is the last record in the file to be pro- 

r 

cessed, the program transfers the contents of the AUX array (which con- 
tains the last record) into work array P. IEXIT is set to 1 and EXT 
called to process the last record. When EXT returns program control to 
MAIN, IEXIT indicates that processing of this file has been completed. 

The program then calls subroutine PRINT to print the irregularity 
table. The analysis of the extracted signal continues by calling the 
statistical analysis subroutine STAT or the spectral analysis subroutine 
EFT. When the analysis is completed for this file of data, control is 
returned to MAIN which loops back to read the instructions for the next 
file to be processed. 

Subroutines 

A description of the subroutines called in the program is given in 
this section. The only subroutines not listed are EOURT and PLOT, since 
they were used in a package furnished by the computer department. A 
Fortran list is included after the description of each subroutine. 

SUBROUTINE RID (13, 14, IC) 

This subroutine reads the identification file which is the first 
file on the tape. This routine uses 3 subroutines that are especially 
written for the University of Alabama IBM-360 Model 50 computer. The 
first, NTRAN, actually reads the tape and stores the data into an array. 
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The data is then transferred from this array and decoded for system 
compatibility using utility subroutines MOVE and TRNSL. This routine 
would probably require modification or rewriting if it-.’were'^used on 
another machine. 


SUBROUTINE REDREC (FILES , IREC, FLICK, IBUF, TIME, N, NOSCAN, NOCHNO, 

ICHNO) 

Subroutine REDREC is called by MAIN to read and reformat the data 
from the magnetic tape. REDREC calls on the special utility subroutine, 
NTRAN to read the tape. Utility subroutines MOVE and TRNSL are called 
to convert the 7 track tape output into the byte system. REDREC must 
unpack from the multiplexed data the desired channel and convert the 
binary code to conventional base ten numbers. It must also keep up with 
the record and file number that it is reading. Provisions were made in 
REDREC to indicate read errors that might occur in NTRAN. This sub- 
routine would probably require modification if used on an other machine. 


Argument Variables 


FILES — Number of files to be processed. 


IREC - Record number as counted by REDREC. 
FLICK - File number as counted by REDREC. 


IBUF - Array containing raw data. 


TIME - Not used. 


N - Number of data points per record (either 200 or 1000) . 

NOSCAN - Number of scans per record per channel (either 200 or 1000). 

NOCHAN - Number of channels multiplexed on the tape. 

Channel number to be processed. 


ICHNO - 
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FORTRAN 

0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 
0019 
0020 - 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 

0051 

0052 

0053 

0054 

0055 


IV G LEVEL 1, BOB 4 RID DATE = 70122 02/35/57 

SUBROUTINE RID (1 3,14 , IC) 

INTEGER IB (5) ,IA (3) , BUF (501) , FLCNT , TTB (3) , BLK, FILES , FILCK , 

*IBUF (1) ,IC (1) 

H=1 

11=0 

12=0 

5 CALL NTRAN ( 1 , 2,3, IA, K, 2,-2004 , BUF, L, 22) 

2 IF (K+ 1) 1,2,3 

3 CALL MOVE (IB,1,IA,1,4) 

CALL TliNSL (IB, 4, TTB) 

DATA TTB/' 012 34 56789’/ 

DATA BLK /' '/ 

IB (2) =BLK 

CALL HOVE (IB(2),1,IA,5,1) 

CALL TRNSL (IB (2) ,1,TTB) 

IB (3) =BLK 

CALL HOVE (IB(3) ,1,IA, 6,1} 

CALL TRNSL (IB (3) , 1, TTB) 

IB (4) =BLK 

CALL MOVE (IB {4} , 1,1 A, 7, 2) 

CALL TRNSL (IB (4) ,2, TTB) 

CALL HOVE (IB (5) , 1,IA,9,4) 

CALL TRNSL (IB (5) ,4, TTB) 

GO TO (6,7) , H 

6 DO 14 1=1,5 

14 IC (I) =IB (I) 

13=11 

14=12 

RETURN 

I TF (K.EQ.-2) GO TO 4 
12=1 

CALL NTRAN (1,22) 

GO TO (6,7) ,M 

4 11=1 

CALL NTRAN (1,22) 

GO TO (6,7) , M 

ENTRY REDSEC (FILES, IREC, FILCK, IBUF, TIME, N, NOS CAN, NOCHAN, ICHNO) 

M=2 

13 IREC=IREC+1 

CALL NTRAN (1,22) 

9 IF (L+1) 8,9,10 

10 TIM E=BUF ( 1 ) 

N= (L~4) / (2+NOCHAN) 

K=3+2*ICHNO 
DO 11 1=1, N 
IBUF (I) =0 
J=0 

CALL MOVE ( J, 4 ,BUF (1) ,K,1) 

CALL HOVE (IBUF (I) ,4, BUF (1) ,K+1, 1) 

IBUF (I) =IBUF (I) + 64* J 

IF (IBUF (I) .GE.1024) IBUF (I) =1 024-IBUF (I) 

II K=K+2*NOCH AN 

CALL NTRAN (1,2, -2004, BUF, L) 

RETURN 

8 IF (L.EQ.-2) GO TO 12 
PRINT 100 
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0056 . 

100 

FORMAT ('0****** READ ERROR ***♦**■) 

0057 


CALL NTRAN (1,22) 

0058 


CALL NTRAN (1,2,-2004, BUF,L) 

0059 


GO TO 13 

0060 

12 

CONTINUE 

0061 


FILCK=FILCK*-1 

0062 


CALL NTRAN (1,22) 

0063 


GO TO 5 

0064 

7 

CONTINUE 

0065 


IREC=Q 

0066 


GO TO 13 

0067 


END 


02/35/57 
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SUBROUTINE EXT (NK, LI, L2, IREG, TOLL, TOL2, T0L3, NCI, IBTYP) 

This subroutine is by far the most complex routine in the program. 
It extracts the intensity amplitudes from the chopped data. A written 
explanation of EXT will not be given due to its complexity. Included 
instead is a flow diagram. It is hoped that the interested reader can 
use the description given in the text along with the flow chart to 
understand the operation of this routine. As an additional aid, the 
important variable names are given a brief description. 



Common Block Variables 

P - Work array containing record of data points being processed. 

AUX - Auxiliary array containing the next record to be processed. 

INDO - Array containing position in the data array from which the 
signal points came. 

JOVF - Number of overflows in HIST. Variable is initialized in EXT. 

JUNF - Number of underflows in HIST. Variable is initialized in EXT. 

Array containing both groups of base points. 


BASE - 
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SIG - Array containing signal points • 

I BLOC “ Array containing position in the data array P from which 

each group of base points came. 

L, M — When EXT calls subroutine LIMIT, it gives it the initial 
position L and the final position M LIMIT is to scan in 
the data array. 

XL90 - This is the result of calling LIMIT and is the criteria 
for selecting base points. 

XL10 - Also the result of LIMIT and is the criteria for selecting 
signal points. 


Labeled Common Variables 

IPRINT - If IPRINT is -other than 0, EXT will print the record if any 
irregularities occur while the record is being processed. 

If IPRINT = 0 it will avoid printing. 

IRTN - Place keeper for subroutine EXT. The subroutine may be in 
any part of its cycle when it completes a record since data 
is continuous from record to record. IRTN is -set to a number 
corresponding to the exit point in the routine when it re- 
turns to the MAIN routine for a new record. When subroutine 
EXT is recalled, a computed GO TO statement keyed to IRTN 
returns control to the phase EXT was previously in. 

IE - Array passed to subroutine print which contains the errors 

accumulated for each record. 

IBCNT - Array containing the number of base points in each group. 

ISCT - Number„of ^signal points. 


Important .Internal Variables 


IRUN - IRUN = 1 indicates routine in start up cycle. Converse for 
IRUN =0. 


INDX - The value of this variable indicates the position (1-200) in 
the record being processed. 

IDINX - -The' number of cpositxons^'subroutirie: EILL 'mustnplace' zeros ^due 
to irregularities which cause data points to be skipped. 

NP - Count of unsuccessful passes through signal and base search 

cycle . 

IBC - Indicates which group of base points are being searched for. 
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LC - Indicates loop condition. LC = 0 means routine is the 

base point search phase; LC = 1 indicates signal point 
search. 

INP - Counts the times the error recycle phase is entered. 

SUBROUTINES HIST (BA, IDUM, NCI, IBTYP , I RUN) 

This routine is called by EXT to compute the amplitude of the 
chopped wave and constructs a histogram with the results. The routine 
is designed to use numbers between 0 - 1000 but can be easily modified 
to handle a larger range. The histogram is stored in a common array 
to be used by the statistical analysis subroutine STAT. 


Argument Variables 


BA - Average of the base points. 


NCI - 
IBTYP - 


IRUN - 


Number of class intervals for histogram. 

Determines the method to be used to calculate the amplitude. 
If IBTYP = 0, amplitudes are calculated by taking the 
difference between the signal points and the average of the 
base points in both groups. If IBTYP = 1, calculation will 
be the difference between the first signal point and the last 
base point of group 1, and the difference between the last 
signal point and the first base point of group 2. If 
IBTYP = 3, calculation is performed only on the last base 
point of group 1 and the first signal point. 

If equal to 1, indicate s' thaf'EXTiis ‘.ini itsf firfetccyfcie. 


Common Block Variables 


NCNT - Array containing frequency for each class interval. 

JOVF - Number of overflows or excessively large numbers resulting 
from classifications of amplitudes. 

JUNF - Number of underflows or very small numbers resulting from 
classification of amplitudes. 

BASE - Array containing two groups of base points. 

SIG - Array containing signal points . 
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FORTRAN 

0001 

0.002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 

0049 

0050 
005.1 
0052 


IV G LEVEL 1, NOD 4 EXT DATE = 70122 

SOBRODTIME EXT (NK,L1 , L2‘, IREC, TOL1 , TOL2 ,TOL3 ,N CI,IBT YP) 
COHBOM P (1000) ,AUX (1000) ,IND0 (10) ,NCNT (300) ,JOVF,JUNF, 
♦ BASE (2,10) ,SIG (10) ,IBLOC (2,10) ,1, H, XL90, XL 10 
CQHNOH/AA A/IPHIHT 
COHHON/RLC/IHTN 
COH BON/PR/IE (7, 300) 

COHHON/CBA/IBCNT (2) ,ISCT 
. IREC=IREC-H 
IST=IST + 1 
IHP = 0 
IHO=1 
INDX= 1 

GO TO (150,11,12,150,776) ,IRTN 
150 IRUM=1 
IRP=0 
IST=0 
L=1 
H=20 

IF (IRTN.EQ.4) 'GO TO 209 
DO 1050 1=1,300 
1050 NCHT(I}=0 
JOVF=0 
JUNF=0 
NDATA=0 

209 CALI II HIT (TOL 1 ,TOL2 , NK) 

C SEARCH FOR FIRST BASE POIHT- HOID VALUE OF IHDX 

776 DO 1 I=L,H 

IF (I.GT.NK) GO TO 2000 
IF (P (I) -XL90) 1,3,3 
3 IHDX=I 
GO TO 4 

I CONTINUE 

IE (1, IREC) =IE (1 , IREC) +1 
IRTH=4 

C RETURN TO DRIVER FOR NEXT RECORD 

RETOEN 
2000 L=1 

M=fi-NK 
IRTN=5 
RETURN 
4 H=INDX+15 
L=IHDX" 

CALL LIBIT (TOL1,TOL2,NK) 

IDIHX=IABS (INDX-INO) 

CALI FILL (IDINX) 

INO=INDX 

IBC=1 

IBCNT (1) =0 
IBCNT (2) =0 
ISCI=0 
10 LC=0 
NP=0 

C BASE POINT SEARCH 

II IF (P (INDX) -XL90) 12,13,13 
13 IBCNT (IBC) =IBCNT (IBC) +1 

I=IBCNT (IBC) 


02/35/57 
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FORTRAN 


0053 

0054 

0055 

0056 

0057 

0058 

0059 
0060- 
0061 
0062 

0063 

0064 

0065 

0066 

0067 

0068 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 
0081 
0082 

0083 

0084 

0085 

0086 

0087 

0088 

0089 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0098 


IV G LEVEL 1, MOD 4 EXT DATE = 70122 02/35/57 

C STORE BASE POINT 

BASE (IBC, X) = P (IHDX) 

IBLOC(IBC,I)=INDX 
IF (IBCNI (IBC) —10) 14,401,401 
401 IBTN=4 

IE{2,IREC}=IE (2,IREC> +1 
IDINX=NK-INDX 
- CALL FILL (IDINX) 

RETURN 

14 INDX=INDX + 1 ' 

INO=INDX 

IF {INDX— HK) 11,11,158 
158 IRTN=2 

C NORMAL RETURN TO DRIVER FOR NEXT RECORD 

RETURN 

C SIGNilL POINT SEARCH 

12 IF (P {INDX) —XL 10) 20,20,15 

15 IF(LC) 16, 16,21 

16 NP= NP+ 1 

C CHECK NUMBER OF UNSUCESSFUL PASSES 

IF (NP-L1) 1944, 1944, 208 

1944 IF(ISCT.EQ.O)' GO TO 14 
GO TO 70'\ 

21 IF (P (INDX) — XL90) 22,30,30 
30 IF (IRUN) 10,10,31 

C SET BASE COUNT INDX CONDITION 

31 IF (IBC-1) 34,33,34 

33 IBC==2 

GO TO 10 

34 IBC=1 

20 IF (IRUN) 23,23,24 
C CHECK LOOP CONDITION 

23 IF(LC) 40,40,29 

24 IF(IBC-I) 25,25,40 

25 1C=1 
NP=0 

29 ISCT=ISCT + 1 
C STORE SIGNAL POINT 

SIG (ISCT) =P (INDX) 

INDO (ISCT) =INDX 
IF(ISCT-IO) 26,501,501 
501 IHTN=4 

IE(3,IREC) =IE (3,IREC) +1 

IDINX=NK-INDX 

CALL FILL (IDINX) 

RETURN 

26 INDX=INDX+1 
INO=INDX 

IF (IHDX— NK) 12,12,162 
162 IHTN=3 

C NORMAL RETURN TO DRIVER FOR NEXT RECORD 

RETURN 

22 NP=NP+1 
If'(NP-LI) 26,26,208 

C CHECK NUMBER OF BASE AND SIGNAL POINTS 

40 IF (IBCHT (1) -LT.TOL3) GO TO 222 


0099 
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FORTRAN IV G LEVEL 1, HOD 4 EXT DATE = 70122 


0100 


IF (IBCNT (2) .LT.TOL3) GO TO 222 

0101 


IF (ISCT.LT.TOL3) GO TO 224 

0102 


GO TO 100 

0103 

222 

IE(4,IREC) =IE (4, IREC) f-1 

0104 


GO TO 555 

0105 

224 

IE (5, IREC) =IE (5, IREC) +1 

0106 


GO TO 555 


C 

ENTER ERROR RECYCLE 

0107 

208 

IE (6 , IREC) =I£ (6, IREC) *1 

0108 

555 

L=INDX 

0109 


M=INDX+20 

0110 


IRUN=1 

0111 


IF (IPRIHT-. EQ-0) GO TO 742 

0112 


PRINT 666, IREC, INDX 

0113 

666 

FORMAT (IX, * IREC= ',110,1 OX, » INDX=‘ ,110) 

0114 


IF (IREC-EQ.IRP) GO TO 742 

0115 


PRINT 333, (P (I) ,I=1,NK) 

0116 

333 

FORMAT ( IX,/, ( IX, 1 0F1 0.3) ) 

0117 


IRP=IREC 

0118 

742 

CONTINUE 

0119 


INP=INP+1 

0120 


IF (IHP-5) 209,921,921 

0121 

921 

IE (7, IREC) =IE(7,IREC) *1 

0122 


IBTN = 4 

0123 


IDINX=NK-INOX 

0124 


CALL FILL (IDINX) 

0125 


RETURN 


C 

POINT SEARCH CYCLE COMPLETE 

0126 

100 

SUM 1=0. 0 

0127 


SUM2=0.0 

0128 


I=IBCNT (1) 


C 

PREPROSSING FOR HISTOGRAM FOLLOWS 

0129 


DO 101 J=1,I 

0130 

101 

S0M1=SUM1+ BASE(1,J) 

0131 


I=IBCNT (2) 

0132 


DO 102 J=1 ,1 

0133 

102 

SUM2=SUM2 +BASE (2,J) 

0134 


BA= (SUH1 +SUM2) / (IBCNT ( 1 ) +IBCNT (2) ) 

0135 


CALL AMPX (IHUN) 

0136 


CALL HIST (BA, NK , NCI , IBTYP , IRUN) 

0137 

43 

IF (IBC— 1) 44,45,44 

0138 

45 

IBC=2 

0139 


GO .TO 46 

0140 

44 

IBC=1 

0141 

46 

ISCT=0 

0142 


IBCNT (IBC) =0 

0143 


LC= 1- 

0144 


NP=0 

0145 


IRUN=0 

0146 


INDX2=INDX+9 

0147 


L-INDX 

0148 


H=INDX2 

0149 


CALL LIMIT (TOL1,TOL2,NK) 

0150 


GO' TO 12 


C LOOK INTO NEXT RECORD FOR SIGNAL POINT 

70 IF {IBCNT (IBC) -L2) 14,72,72 


02/35/S7 
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FOBTRAN 

0152 

0153 

0154 

0155 

0156 


IV G LEVEL 1, HOD 4 EXT DATS = 70122 

72 L=IHDX 

H— IHDX + 9 

CALL LIMIT (TOX 1 ,TOL2 , UK) 

IF (P (IHDX) —XL 10) 20,20,14 
END 


02/35/57 
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FORTRAN 

0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 


IV G LEVEL 1, MOD 4 HIST DATE = 70122 02/35/57 

SUBROUTINE HIST (BA ,IDUM , NCI , IBTTP ,IRUN) 

COMMON P(1000) , AUX (1000) ,IND0(10) ,NCNT(300) .JOVF.JUNF, 

*BASE{2, 10) ,SIG (10) ,IBLOC (2,10) , KK, HM, XL90 , XLl 0 
DIMENSION Y (40) , X (40) , LOC (4 0) , A(15), V(15), B (200) 
COMMON/CBA/IBCNT (2) ,N 
C COMBINE BASE 1 AND BASE 2 

IF (IBTTP) 100, 100,200 
200 IE (IRUN.EQ. 1) 4mHH=1 
KKK=1 
LLL=1 

IF (MMM.EQ. 1) LLL=2 
IF(MMM.EQ.O) KKK=2 
100 IJ=N 

DO 1 1=1, N 

IF (IBTYP.EQ.O) GO TO 202 

400 IF(I-I) 402,401,402 

401 NTEKP=IBCNT (KKK) 

BA=8ASE (KKK, NTEMP) 

GO TO 202 

402 IF(I-IJ) 1,403,1 

403 BA=BA$E (LLL, 1 ) 

202 J=BA-SIG(I) 

J= (NCI*J)/1000 + 1 

IF (J-1) 2,3,3 

2 JUNF=JUNF+1 
GO TO 1 

3 IF (3-300) 5,5,4 

4 JOVF=JOVF+1 
GO TO 1 

5 NDATA=NDATA+ 1 
NCNT(J)=NCHT (J) +1 

IF (IBTTP. EQ. 3) GO TO 50 
1 CONTINUE 

IF (BHM.EQ .0) MHH=1 
IF (MMM.EQ. 1) MMH=0 
50 RETURN 
END 
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SUBROUTINE STAT (NTAPE , NCH, NFILE, NCI, I FLAG) 

STAT performs the statistical analysis by using the intensity 
histogram constructed by HIST. A log-amplitude histogram is generated 
from the intensity histogram to perform log-normal tests. The mean, 
standard deviation, skewness and kurtosis are calculated for both the 
intensity and log-amplitude data. In addition, the cumulative probability 
is calculated and a chi-square test made on both sets of data. 



NCI - Number of class intervals. 

IFLAG - Indicates type statistical calculations. 


Common Block Variables 

NCNT - Array contains histogram. 

JOF — Number of overflows. 

JUF — Number of underflows. 


SUBROUTINE LIMIT (T0L1, T0L2, NK) 

This routine is called by subroutine EXT to calculate the criteria 
for determining if a data point is a base point or a signal point or 
neither . LIMIT has the capability of looking into the next record if it 
is called near the end of the record being processed. 




FORTRAN IV G LEVEL 1, MOD 4 


STAT 


DATE = 70122 


02/35/57 


0001 



SUBROUTINE STAT (NTAPE , NCH, N FILE, NCI , IF LAG) 


C 


IFLAG = I NO CHI SQUARE TEST 


C 


2 CHI SQUARE TEST ON NORMAL DISTRIBUTION 


C 


3 CHI SQUARE TEST ON LOG NORMLL DISTRIBUTION 


c 


4 CHI SQUARE TEST ON BOTH 

0002 



COMMON P(1000) ,AUX (1000) ,INDO(10) ,!ICNT1 (300) , JOF , JDF,BASE (2,10) 



4 

“PIG (10) , IBLQC (2,10) , KK, MM , XL90 , XL 10 

0003 



COHHON/UUU/NPNCH 

0004 



DIMENSION Y (300) ,XLN (300) ,Q (4) ,R (4) 

0005 



DIMENSION NCNT (300) 

0006 



DO 60 1=1,300 

0007 


60 

NCNT (I) =NCNT1 (I) 

0008 



C=1000/HCI 


c 


FIND THE HIGHEST AND LOWEST CLASS INTERVAL 

0009 



DO 1 1=1, NCI 

0010 



IF (NCNT (I)) 1,1,2 

0011 


2 

ILO=I 

0012 



GO TO 3 . 

0013 


1 

CONTINUE 

0014 


3 

DO 4 1= ILO , NCI 

0015 



IF (NCNT (I) ) 4,4,5 

0016 


5 

IHI=I 

0017 


4 

CONTINUE-. 


c 


FIND NUMBER OF POINTS AND FLOAT NCNT 

0018 



N=0 

0019 



DO 6 1= ILO, IHI 

0020 



Y (I)=NCNT (I) 

0021 


6 

U=H+NCNT (I) 


c 


PRINT HEADINGS 

0022 



WRITE (6,101) NTAPE, NCH, NFILE 

0023 



WRITE (6,1 02) 


c 


DEFALT DUE TO TOO FEW CLASS INTERVALS 

0024 



NN= IHI-ILO 

0025 



IF (NN-10) 50,50,51 

0026 


50 

WRITE (6,110) 

0027 



RETURN 


c 


FIND AVERAGE AMPLITUDE 

0028 


51 

XN=N 

0029 



AVE=0.00 

0030 



DO 8 1= ILO, IHI 

0031 



XI=I 

0032 


8 

AVE=AVE+Y (I) * (XI-0.5) *C 

0033 



AVE=AVF./XN 


c 


COMPUTE CUMULATIVE PHOB ABILITIEC AND LOG AMPLITUDES 

003*1 



SUM=0. 0 0 

0035 



DO 10 I=ILO,IHI 

0036 



xr=i 

0037 



xi= (xi-o.S)*c 

0038 



XLN (1) =0. 5+ALOG (XI/AVE) 

0039 



SUH=SUM+Y (I) 

0040 



CP=SUB/XN 

0041 


10 

W RITE (6,1 03) XI, XLN (I) ,Y (I) ,CP 

0042 



WRITE (6, 111) JOF, JUF 


c 


COMPUTE MOMENTS ABOUT THE MEAN 

0043 



XLA=0. 00 

0044 



DO 20 I=ILO,IHI 
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0045 

20 

XLA=XLA+Y (I) *XLN (I) 

0046 


XLA=XLA/XH * 

0047 


DO 21 1=2,4 

0048 


Q(I) =0.00 

0049 

21 

R (I)=0.00 

0050 


DO 22 I=ILO, IHI 

0051 


XI=I 

0052 


DO 22 J=2, 4 

0053 


Q (J) =Q (J) +Y (X) * ( (XI-0. 5) *C- AVE) 

0054 

22 

R (J) =R (J) + Y (T) * (XLN (I) -XLA) 

0055 


DO 23 J=2, 4 

0056 


Q(J)=Q(J)/XN 

0057 

23 

R ( J) =R (J) /XN 

0058 


NT=IHI— ILO+ 1 

0059 


WRITE (6, 104) NT 

0060 


SIG=SQRT (Q (2) ) 

0061 


SIGL=SQRT (R {2) ) 

0062 


SKEW=0.5*Q (3)/(SIG**3) 

0063 


SKL = 0.5*R (3)/ (SIGL**3) 

0064 


XRUR= ( (Q (4) / (Q (2) **2 ) ) -3 . 0) /2 . 0 

0065 


XKURL = ((R(4)/(R(2)**2))-3,0)/2,0 


c 

PRINT HOMENTS 

0066 


WRITE (6,105) AVE, SIG, SKEW, XKUF, XLA, SIGL,SKL, XKURL, H 

0067 


IF (NPNCH .EQ. 0) GO TO 810 

0068 


PUNCH 800, NT APE, NCH , H FILE, A V£ , SIG ,SIGL, XLA 

0069 

800 

FORHAT (A4,I2,I4,4E14.4) 

0070 

810 

CONTINUE 

0071 


GO TO (31,32,42,32) , IFLAG 

0072 

31 

RETURN 


c 

Q 

NO CHI SQUARE TEST REQUESTED 


c 

CHI SQUARE TEST 

0073 

32 

CALL CHI (CSQ, Y,ILO, IHI, C, NUSE, AVE, XLA, SIG ,XN,Q,SIG) 

0074 


XN=CSQ 

0075 


WRITE (6,106) CSQ , NUSE 


C 

PRINT CHI SQUARE NORHAL 

0076 


GO TO (31,31,42,42) , IFLAG 

0077 

42 

CALI CHI(CSQ,Y,ILO, IHI, C, NUSE, AVE, XLA, SIG ,XN,1,SIGL) 

0078 


AXLN=CSQ 

0079 


WRITE (6,106) CSQ, NUSE 


C 

PRINT CHI SQUARE LOG NORHAL 

0080 


PUNCH 999,NTAPE, NCH, N FILE, SI GL, SKEW ,5 KL,XN , AX LN 

0081 

999 

FORHAT (A4, 1 1,12, 5E 13. 6) 

0082 


RETURN 

0083 

101 

FORHAT (' 1 ' ,5X, ‘TAPE NUHBER A4 , 5X, ' TRACK* ,13 ,5X, » FILE* ,13) 

0084 

102 

FORM AT (’ O' ,16X,' AMPLITUDE' ,1 OX, 'LOG AMPLITUDE* ,12X, 'COUNT' , 
1 8X, ' CUMULATIVE PROBABILITY'/) 

0085 

103 

FORHAT (7X,4E21. 6) . 

0086 

104 

FORHAT(«0',10X, 'NUMBER OF CLASS INTERVALS = ',I6) 

0087 

105 

FORHAT ('O' , 17X, ' AVERAGE* ,9X, ' STANDARD DEVIATION * , 8X SKEW NESS • , 


*1 3X , 'KUETOSIS '/7X, 4E2 1 . 6/7X ,4E21.6/'0',10X,* NUMBER OF DATA POINTS 



*,110) 

0088 

106 

FORHAT ('O' ,10X,' CHI SQUARE= ' , El 4 . 6/1 1 X ' NUHBER OF CLASS INTERVALS 
1USED =• , 15) 

0089 

110 

FORMAT {'0*,5X,' TOO FEW CLASS INTERVALS '/ 

1 ' 0 • , 5X , * EXECUTION OF STATISTICS CALCULATION SUSPENDED’) 
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FORTRAN IV G LEVEL 1, HOD 4 


STAT 


DATE = 70122 


0090 


0091 


111 FORMAT('0*,5X,'NUHBER OF OVERFLOWS' ,16/ 
*6X, 1 NUMBER OF U NDERFLOWS ' ,15) 

END 


02/35/57 
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Common Block Variables 


Gives the position (value of : INDX) in the record at the 
time LIMIT is called. 

This variable is the sum of KK and the number of points 
LIMIT is to scan. 


XL90 - 


The resulting criteria for base point selection. 


XL10 - 


The resulting criteria for signal point selection. 


SUBROUTINE AMPX (IRUN) 


This routine takes the signal and base points extracted by EXT and 
computes the amplitude of the square wave for the spectral analysis. The 
amplitudes are calculated by taking the difference between the last base 
point in the. first group and the first signal point, and the difference 
between the last signal point and the first base point in the second group. 
This produces two signal points per group. The points are stored in 
array AMP for use by the spectral analysis routines . 


Argument Variables 

IRUN - Indicates if EXT in startup cycle. 

Common Block Variables 

BASE - Array containing both groups of base points. 

SIG - Array containing signal points. 


Labeled Common Variables 

IBCNT - Array containing number of base points in each group. 

N - Number of signal points 

SUBROUTINE CHI (CSQ, Yl, ILO, IHI, C, NUSE, AVE, XL A, SD, XN, NTYP, SX) 
This routine is called by the statistics subroutine to perform the 
chi-square test for normal and log-normal distributions. 
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FORTRAN 

0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 
0 039 

0040 

0041 


IV G LEVEL 1, MOD 4 LIMIT DATE = 70122 02/35/57 

SUBROUTINE LIMIT (TOX1 ,TOL2 , NK) 

COMMON P{1000) ,AUX(1000) ,INDO(10) ,NCNT{300) ,JOVF,JDNF, 

♦BASE (2, 10) ,SIG{10) ,IBLOC (2,10) ,KK, MM, XL90, XL1 0 
IK=KK 
IM=HH 
ID=0 

IF (IK-NK) 502,501,666 
666 PRINT 3000 

3000 FORMAT (' IK GREATER THAN NK ') 

501 AMAX=P(IK) 

AMI N=AM AX 
ID= IM-NK 
GO TO 381 

502 IF(IH-NK) 360,360,361 
361 ID=IM-NK 

IH=NK 

360 AMAX=P(IK) 

AMIN=AMAX 

DO 350 J=IK,IM 

IF (AMAX— P (J) ) 301,302,302 

301 AMAX=P(J) 

302 IF (AMIN— P (J) ) 350,303,303 

303 A«IN=P(J), 

350 CONTINUE 

IF (ID) 380,380,381 
381 AM AXP=AUX ( 1) 

AMIHP=AHAXP 
DO 650 J=1,ID 

IF (AHAXP— AUX (J) ) 601,602,602 

601 AMAXP^AUX (J) 

602 IF (AMINP-AUX (J) ) 650,603,603 

603 AMINP=AUX(J) 

650 CONTINUE 

IF (AMAX— AMAXP) 700,701,701 

700 AHAX=AMAXP 

701 IF{AMIH-AMINP) 380,380,703 
703 AMIN=AMINP 

380 A=AMAX- AMIN 

XL9 0=AHAX— TOI 1*A 
XL10=AMIN+ TOL2*A 
RETURN 
END 
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FORTRAN IV G 1EYEL 1, MOD 4 ABPX DATE = 70122 


0001 


SUBROUTINE AMPX(IRUN) 

0002 


COMMON P{1000) f AUX (1000) ,INDO(10) ,HCNT{300) ,OOVF,JUNF, 
*B AS E (2 , 1 0) ,SIG(10) ,IBiOC (2 , 1 0) ,KK , MM, XL90 , XL1 0 

0003 


COMMON/CBA/IBCNT (2) ,N 

0004 


COMMON/AP/AHP (10000) , NDATA 

0005 


IF(IRUN.EQ. 1) MHM=1 

0006 


KKK=1 

0007 


UL=1 

0008 


IF(MMH.EQ.I) 1LL=2 

0009 


IF {MHM. EQ- 0) KKK~2 

0010 


IJ=N 

0011 


DO 1 1=1, N 

0012 

400 

IF{I-1) 402,401,402 

0013 

401 

NTEMP=IBCNT (KKK) 

0014 


BA=BASE (KKK, NTEMP) 

0015 


GO TO 502 

0016 

402 

IF(I-IJ) 1,403,1 „ 

0017 

403 

BA=BASE (ILL, 1) 

0018 

502 

NDATA=NDATA+1 

0019 


AMP (NDATA) =BA-SIG (I) 

0020 

1 

CONTINUE 

0021 


IF (MMH.EQ. 0) MHH= 1 

0022 


IF (MMH.EQ. 1) MHH=0 

0023 


RETURN 

0024 


END 


02/35/57 
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FORTRAN IV G LEVEL 1, HOD 4 CHI 


DATE = 70122 02/35/57 


0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 

0031 

0032 

0033 

0034 

0035 

0036 

0037 

0038 

0039 


SUBROUTINE CHI(CSQ,Y1 ,IL0 ,IHI,C, NUSE, A VE,XLA, SD,XH, NTI P, SX) 
DIHENSION Y (300,3) , 71 (300) 

DO 1 1=1,300 

1 Y (1,1) =0.00 
NUSE=0 

KH= (ILO+IHI) /2 
J=1 

Y (J,2) =AVE— 1 0. 0*SD 

GROOP CLASS INTERVALS ON LOW END 
DO 2 1= ILO.Kfl 

Y (J,1)=Y1 (I) + Y (J,1) 

IF (Y (0,1) —5.0) 2,2,3 

3 Y(J,3)=C*I 

NUSE = NUSE + 1 
J = 0 + 1 

Y ( J, 2) = C*I 

2 CONTINUE 

C GROOP CLASS INTERVALS ON HIGH SIDE 

I = IHI 

Y(J,3) = AVE * 10.0*SD 
6 IF(I-KM) 10,10,4 

4 Y (0,1) = Y (J,1) + Y1 (I) 

IF(Y(J,1) r 5.0) 11,11,5 

5 Y (J,2) = C*(I-1) 

NUSE = NOSE + 1 
J = J t 1 

Y (0,3) = C*(I-1) 

111=1-1 

GO TO 6 

C COMPUTE THEORITICAL PROBABILITY 

10 CSQ = 0.00 

DO 30 1=1, NUSE 
XLL = I (1,2) 

XUL=Y (1,3) 

24 CALL SIMP { FTH, XLL ,XUL,2 1 , HTYP, A V£, SX ,XLA) 

C COMPUTE CHI SQUARE 

IF (FTH) 31,31,30 
31 WRITE (6, 100) I 

100 FORMAT ( 1 0 ZERO VALUE OF THEORITICAL PROBABILITY IN', 15,* TH 
1 INTERVAL'/ 6X, 'EXECUTION OF CHI SQUARE TEST DISCONTINUED’) 
RETURN 

30 CSQ=CSQ+((Y (I,1)-XN*FTH)**2)/(XN*FTB) 

RETURN 

END 
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Argument Variables 


CSQ - 

Result of chi-square test. 


Y1 - 

Array containing histogram. 


ILO - 

Lowest class interval in histogram. 


IHI - 

Highest class interval in histogram. 


C - 

Width of class mark. 


NUSE - 

Number of class intervals used by the 

chi-square routine. 

AVE - 

Mean value of amplitudes. 


XLA - 

Mean value of log-amplitudes. 


SD - 

Standard deviation of amplitudes . 


XN - 

Number of data points. 


NTYP - 

Determines if chi-square test will be 
normal test or both. 

run for normal or log- 

SX - 

Log standard deviation. 


SUBROUTINE 

FFT. (DT, FMAX) 



This routine is called by the main program to coordinate the per- 
formance of the spectral analysis. Subroutines FOURT and PLOT are called 
to perform the Fourier transform and to plot the results. FFT will have 
PLOT plot directly from the calculated spectral data array or it will have 
it plot the average of a designated number of points in the array. This 
feature was incorporated to smooth out random variations . 

Argument Variables 

DT - Time between data samples. 

FMAX - Maximum frequency to be used. 

Labeled Common Variables 

AMP - Array containing the time domain signal. This array is 

passed from AMPX. 
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FORTRAN IV G IE VEL 1, HOD 4 FFT 


DATE. = 70122 


0001 


5UBR00TINE FFT(N,DT,FMAX) 

0002 


DIMENSION WORK (2500) 

0003 


COMHON/AP/AMP (10000) , NDATA 

0004 


COMMON/GO/H1.N3,N4,FH,IDO 

0005 


COMMON/NPTS8/N8 

0006 


DIMENSION NN{2) 

0007 


PRINT 555, NDATA 

0008 


. IF(HDATA-8192) 25,56,56 

0000 

25 

NDIFF=8 192~ NDATA 

0010 


CALL FILL (NDIFF) 

0011 

56 

PRINT 81, N8 

0012 

81 

FOR HAT (10X,* NUMBER OF ZEROS ADDED = *,120//) 

0013 


NDATA=8 192 

0014 

555 

FORMAT ( 10X , * NDHBER OF DATA POINTS EXTRACTED =*,120,///) 

0015 


N=NDATA 

0016 

99 

FORMAT (/,' 1FHEQ (HZ) MAGNITUDE ' ) 

0017 


NN ( 1) =N 

0018 


DF=1.0/(N*DT) 

0019 


N2= (FHAX/DF) *2 

0020 


IF(N2.GT.N) N2=N 

0021 


CALL FOURT (AMP, NN,1, -1,0, WORK, 2500) 

0022 


DO 5 J=2, N,2 

0023 


X=AMP (J-1) 

0024 


I=AMP (0) 

0025 

5 

AMP (J- 1) =SQRT (X*X+Y*Y) /N 

0026 


WRITE (6,99) 

0027 


IF (IDO. EQ. 1) GO TO 50 

0028 


CALL PLOT (AMP, DF,DF,N1,N2,N3) 

0029 

50 

F=DF 

0030 


IF (IDO. EQ. 3) GO TO 51 

0031 


PRINT 100 

0032 

100 

FORMAT (/,' FREQ (HZ) MAGNITUDE*) 

0033 


SUH=0.0 

0034 


J=1 

0035 


NC=0 

0036 


DO 10 I=N1,N2,N3 

0037 


HC=NC+1 

0038 


SUM=Af1P(I) +SUM 

0039 


IF (NC.LT.N4) GO TO 10 

0040 


NC=0 

0041 


AMP (U) = SUM/N4 

0042 


SUH=0. 0 

0043 


IF(F.GT.FH) GO TO 11 

0044 


J=J+1 

0045 

10 

F=F+ DF 

0046 

11 

DFP=N4*DF 

0047, 


SF=(DF + DFP)/2.0 

0048 


CALL PLOT (AMP, SF,DFP, 1,J, 1) . 

0049 

51 

CONTINUE 

0050 


RETURN 

0051 


END 


02/35/57 
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NDATA - Number of data points in AMP array. 

— Designates the first point to be plotted from the spectral 

data array by PLOT. 

N3 _ Directs PLOT to skip N3 points between each plotted point 

in spectral data array, 

N4 - Number of points to be averaged when using the averaging 

feature of this routine. 

IDO - If IDO = 1 the 11 average" feature is to be used. If IDO = 3 

the "average" feature is not to be used. 

N8 - Number of points (zeros) added by FILL. 

SUBROUTINE PRINT (NNN) 

This subroutine accepts the error cumulation array from EXT after 

a run is completed. It prints out the error table and other error data. 

Labeled Common Variables 

IE - Array- containing the sum of 7 types of errors for each of 

the 300 records. 

IRCEND - The -number of the last record to be processed. 

SUBROUTINE SIMP (SUM, FLL, FUL, N, NTYP , A, B, C) 

This isra' Simpsons rule integration routine called by CHI. 



FUL - Upper limit. 

N - Number of points - 21. 

NTYP - Determines if routine will compute for a normal test, log- 
normal test or both. 

A - Mean value of amplitudes 

B - Log standard deviation. 

G - Mean value of log-amplitude. 
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FORTRAN 

0001 

0002 

0003 

0004 

0005 

0006 


0007 

0008 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 
0021 
0022 

0023 

0024 

0025 

0026 

0027 

0028 

0029 

0030 


IV G LEVEL 1, HOD 4 PRINT DATE = 70122 02/35/57 

5GBR0UTIHE PRINT (NNN) 

COB HON/PR/IE (7, 30 0) 

COHHON/RAT/IECEND 
DIMENSION IESUM (7) 

PRINT 3 

3 FORMAT (1H1 ,501. » DATA PROCESSING IRREGOLARTIES ' ,/,10X r ' ERROR CODES 
• *FOLLOW ',////, 

*10X,'H0 BASE POINTS FOUND IN BASE SEARCH ********** 1 •/, 

*10X, 'NUMBER OF BASE POINTS EXCEEDS 10 ************ 2 '/, 

*1 OX ,' HUMBER OF SIGNAL POINTS EXCEEDS 10 *********** 2 '/, 

*10X, 'NUMBER OF BASE POINTS INSUFFICIENT *********** 4'/, 

*10X, 'NUMBER OF SIGNAL POINTS INSUFFICIENT ********* 5'/, 

*10X, 'NUMBER OF PASSES EXCEEDS LIMIT LI ************ 6'/, 

*10X, 'NUMBER OF ERRORS IN RECORD EXCEEDS 5 ********* 71) 

DO 555 1=1,7 
555 IESUM (I}=0 
PRINT 1 

1 FORMAT 1 10X , ' ERROR 1 ,11X,'1',13X,«2 , ,1 3X, *3',13X,'4',13X, *5',13X, 

*'6',13X,'7') 

PRINT 100 

100 FORMAT (10X, ' RECORD 1 ,///) 

DO 55 I=1,IRCEND 
DO 55 K= 1 ,7 

55 IESUM {K) =IESUM (K) +IE (K,I) 

DO 2 1= 1 , IRCEND 

DO 12 K=1,7 

IF (IE {K, I) ) 15,12,15 

12 CONTINUE 
GO TO 2 

15 HRITE (6,5) I, (IE (K,I) ,K=1,7) 

5 FORMAT (10X,I3,10X,I4, 10X,I4 , 1 0X, 1 4, 1 0X ,14 , 1 0X, 14 , 10X, 14 , 1 0X , 14) 

2 CONTINUE 
PRINT 20 

20 FORMAT (10X,//,10X,' ERROR CODE' , 1 OX, • NUMBER OF ERRORS') 

DO 50 K=1,7 

50 PRINT2 1 , K , IESUM { K) 

21 FORMAT{10X,I6,14X,I9) 

RETURN 

END 



FORTRAN IV G LEVEL 1, MOD 4- 


SIMP 


DATE = 70122 


02/35/57 


0001 . 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 


0011 

0012 

0013 

0014 

0015 

0016 

0017 

0018 

0019 

0020 

0021 

0022 

0023 

0024 

0025 

0026 
0027 
0029 

0029 

0030 

0031 

0032 

0033 


SUBROUTINE SIMP (SUM ,FLL, FUL, N, NTYP, A , B , C) 

C INTEGRAND FUNCTION REMOVE HHEHCHANGING FUNCTION 

PBF(X,A,S) = (1,'0/(S*5QRT (6-28318) ))*EXP (-0-5* ( ( (X-A) /S) **2)} 
FNP=N-1 

DELX= (FUL-FLL)/FNP 

SUM=0. 0 

SUM 1=0. 0 

SUM2=0.0 

DO 1 1=1, N 

FK=I-1 

X=FK*DELX+FLL 

C CALL FOR INTERGRAND SOBROUTIHE HERE 

C CALL FOR INTERGRAND SUBROUTINE HERE 

IF(NTYP) 101,101,110 

101 VAL=PBF (X, A,B) 

GO TO 102 

110 IF(X) 20', 20, 100 
20 VAL=0.00 
GO TO 102 

100 XX=0.5*ALOG (X/A) 

VAL=PBF(XX,C,B) 

VAL=0.5*VAL/X 

102 CONTINUE’. 

C 

IF(I.EQ.I.OR.I.EQ.N) GO TO 2 
J=M0D(I,2) 

IF (J) 3, 4,3 

3 SUH1=SUM1+VAL 
GO TO 1 

4 SUS2=SUM2+VAL 
GO TO -1 

2 SUM=SUM+VAL 
1 CONTINUE 

SUM=SUB+2. 0*SUM1 * 4 . 0*SUM2 

SUM=SUM*DELX/3.0 

RETURN 

END 
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SUBROUTINE FILL (I) 

This routine is called by subroutine EXT in cases where data points 
are discarded due to irregularities. FILL places zeros into the omitted 
positions. 

Argument Variables 

I - Number of data points discarded. 

Labeled Common Variables 

AMP - . Array containing extracted data. 

NDATA - Number of data points in AMP. 

N8 - Number of zeros added by FILL. 



FORTRAN 

0001 

0002 

0003 

0004 

0005 

0006 

0007 

0008 

0009 

0010 


IV G LEVEL 1, HOD 4 FILL DATE = 70122 

SUBROUTINE FILL (I) 

COHBON/AP/AHP (10000) , NDATA 

COMHON/NPTS8/N8 

J= (1+3) /6 

N8=N8+J 

DO 1 K=1,J 

ND AT A=ND ATA + 1 

AHP (NDATA) =0, 

RETURN 

END 


02/35/57 



ERROR COOES FOLLOW 


DATA PROCESSING I RREGUL ART1 ES 


NO BASE POINTS FOOND IN BASE SEARCH ********** 1 

NUMBER OF BASE POINTS EXCEEDS 10 ************ 2 

NUMBER OF SIGNAL POINTS EXCEEDS 10 *********** 3 

NUMBER. OF BASE POINTS INSUFFICIENT *********** 4 

NUMBER OF SIGNAL POINTS INSUFFICIENT ********* 5 

NUMBER OF PASSES EXCEEDS LIMIT LI ************ 6 

NUMBER OF ERRORS IN RECORD EXCEEDS 5 ********* 7 

ERROR '1 2 3 

RECORD 


4 


5 


6 


2 

232 

240 

241 
25 2 
253 

264 

265 

276 

277 
238 
289 
300 



0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 


ERROR COOE 
1 
2 

3 

4 

5 

6 
7 


NUMBER OF ERRORS 
0 
1 
0 

12 

0 

0 

0 


VO 


ERROR TABLE 


OOOOOOOOOOOOO 



TAPE NUMBER 3397 


TRACK 4 


FILE 2 


AMPLITUDE 

LOG AMPLITUDE 

COUNT 


CUMULATIVE PROBABILITY 

0.235000E 

03 

-0.347886E 00 

0 . LOOOOOE 

01 

0. I20077E— 03 

0.245000E 

03 

-0.327O5OE 00 

0.200000E 

01 

0 i 360 230 E— 03 

0 • 255000E 

03 

— 0.307047E 00 

0.20Q000E 

01 

0.600384E— 03 

0. 26 5 000 E 

03 

— 0 . 287B14E 00 

0.0 


0.600384E— 03 

0.2750Q0E 

03 

-0.269293E 00 

0.0 


0.600384E— 03 

0.285000E 

03 

-0.2514346 00 

0.200000E 

01 

0 . 840538E— 03 

0 . 295000 E 

03 

—0 .2341916 00 

0 .100000E 

01 

0.960615E— 03 

0.305000E 

03 

-0 .2175236 00 

0.300000E 

01 

0.132085E— 02 

0.315000E 

03 

-0.201393E 00 

0 .6 OOOOOE 

01 

0.204 13 IE— 02 

0.325000E 

03 

— 0.185766E 00 

0.800000E 

01 

0.300I92E-02 . 

0 • 335000E 

03 

— 0 . 170614E 00 

0.400000E 

01 

0.348223E— 02 

O.345O0OE 

03 

-0.155907E 00 

0.170000E 

02 

0.552353E— 02 

0.355000E 

03 

— 0. 141620E. 00 

0 .190000E 

02 

0. 7804996-02 

0.365000E 

03 * 

— 0 • 127730E 00 

0.210000E 

02 

0 • 103266E— 01 

0.375000E 

03 

-0.114216E' 00 

0 .3900006 

02 

0.150096 E-01 

0. 385000E 

03 

-0 .1010576 00 

0.660000E 

02 

0 • 229347E— 01 

Q.395Q00E 

03 

-0.882359E-01 

0.7200006 

02 

0.315802E— 01 

' 0.405000E 

03 

-0.757353E-01 

0.127000E 

03 

0.468300E— 01 

0.4I5000E 

03 

-0.635396E— 01 

0.178000E 

03 

0. 6820366—01 

0.425000E 

03 

-0.516343E-01 

0.286000E 

03 

0 . 102546E 00 

0.435000E 

03 

— 0.40005 8 E— 01 

0.381000E 

03 

0.148295E 00 

0.4450006 

03 

-0.286417E— 01 

0.513000E 

03 

0.209894E 00 

0.455000E 

03 

-0 .175301E-01 

0.739000E 

03 

0.298631E 00 

0.465000E 

03 

-0.6660 13E— 02 

0.929000E 

03 

0.410182E 00 

0.475000E 

03 

0.397811E— 02 

0.1 10800E 

04 

0. 543228E 00 

0.485000E 

03 

0.143953E— 01 

0. 122700E 

04 

0.690562E 00 

0.495000E 

03 

0. 245999 E-OL 

0 .1 17400E 

04 

0. 831532E 00 

0.505000E 

03 

0.346000E— 01 

0.855000E 

03 

0.934198E 00 

0.515000E 

03 

0.4440 44 E-01 

0. 434000E 

03 

0.986311E 00 

0.525000E 

03 

0.540202E— 01 

0.105000E 

03 

0. 998919E 00 

0.535000E 

03 

0.634542E-01 

0.900000E 

01 

O.IOOGOOE 01 

NUMBER OF OVERFLOWS 

0 





NUMBER OF UNDERFLOWS 0 





NUMBER OF CLASS INTERVALS = „ 31 




AVERAGE 


STANDARD DEVIATION 

SKEWNESS 

KURTOSIS 

0.471235E 

03 

0.323860E 02 

-0.6461816 

00 

0. 163909E 01 

-0. 127374E- 

-02 

0.364458 E-01 

— 0-903620E 

00 

0.35 3928E 01 

N UMBER OF DATA 

POINTS 

8328 




CHI SQUARE® 0. 

■ 542616E 

06 




NUMBER OF CLASS INTERVALS USED = 23 




, 



TYPICAL COMPUTER OUTPUT FOR 

CHI SQUARE- 0. 

534447E 

07 

STATISTICAL 

ANALYSIS 

NUMBER OF CLASS INTERVALS USED = 23 






TAPE NUMBER=3397 IDCAL=3 SASE=1 NUMBER OF CHANN E LS=0 5 

TYPE BASE CALCULATION a**##*#*#***^##**#***#*:**********#***# 
TYPE STATISTICS CALLED FOR fc**************!******** *$*##*#*** 
NUMBER OF CLASS INTERVALS mw^tfT##***’**#*)!^* 1 !'*#*********:!:*** 
NUMBER. OF SCANS PER RECORD if,*#******#**^*#*#***#*:*#'!"****#*** 
NUMBER OF CHANNELS ON TAPE ****#$4=fc*#*##*=**>i!#*#Y***>i‘**:4<‘4##* 
TOLERANCE ON BASE LIMITS 

TOLERANCE ON SIGNAL LIMITS ***#-*$ ***#$**#*:****3***«*s#**#* #* 
NUMBER OF SIGNAL AND BASE .POINTS REQUIRED PER CYCLE ******** 
NUMBER OF PASSES REQUIRED TO ABORT CYCLE ******************* 
NUMBER OF', BASE POINTS REQUIRED FOR SIGNAL POINT SEARCH ***** 
PUNCHED OUTPUT FLAG **************************************** 
PR"[NT ERROR TABLE FLAG ************************************* 
TYPE ANALYSIS ROUTINE rfANTEO ******************************* 
INFORMATION FOR FFT ROUTINE ******************************** 


TYPICAL PARAMETER VALUES 


NUMBER OF SCAN S=0200 


* **** 

IBTYP 

***** 

* 

2 


I FLAG 

***** 


4 

***** 

NCI 

***** 


100 

***** 

NOSCAN 

***** 


200 

***** 

NOCHAN 

***** 


5 

***** 

TOLL 

***** 


0.050 

***** 

TOL2 

***** 


0. 100 

***** 

T0L3 

***** 


3.000 

***** 

U 

***** 


10 

***** 

L2 

***** 


5 

.***** 

NPNCH 

***** 


1 

***** 

I PRINT 

***** 


0 

***** 

I SS8 

***** 


2 

***** 

FMAX 

***** 


40.000 

***** 

DT 

***** 


0.006 

***** 

FM 

***** 


30.000 

***** 

N1 

***** 


3 

***** 

N 3 

***** 


2 

***** 

N4 

***** 


5 

***** 

IDO 

***** 


5 
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APPENDIX B 

This program generates a set of N random numbers having a log- 
normal distribution and a pre-selected mean and standard deviation. The 
program is in the form of a FORTRAN IV subroutine. 

Theory : By definition a log-normal random deviate is one whose logarithms 

are normal random deviates. Thus if (X^ is a set of log-normal random 
numbers then there must exist a set of normal random numbers (y.) related 

i 

to the X_^ by 

y i = In X ± B1 

Equation B1 may be generalized by the addition of appropriate scaling 
factors; i.e., we may let 

y^ = a In X_^ 4- b ■ B2 

Now by choosing the mean and variance of the (y^) and the values of the 

scale factors a and b it is possible to generate a set of (X_^) having any 

desired mean and variance from a set of normal deviates (y ) . Solving 

B2 for X. we have 
i 

y. - b 

, X. = exp (— ) B3 

X 3 . 

Since we wish to specify only two parameters, viz., the mean and standard 
deviation of the (X_^) it seems reasonable to assume that we will need 
only two parameters in equation B3. We therefore let a = 1 and take the 
mean of the (y^) to be zero. B3 then becomes 

X. = exp (~b) exp (y^ 


B4 
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taking the average of both sides of equation. B4 we have 

X = exp (-b) exp (y ± ) B5 

and also taking the second moment of (X_p about zero 

~2 

X = exp (-b) exp (2y ) 336 

the averages of the exponential functions in equation B5 and B6 can 
be evaluated easily 


2 - 1/2 

exp(ny ± ) = (2irt ) 


f X 2 

exp(ny)* exp(y /2o)dy 

. 

-x 


B7 


Combining equation B5, B6 and B7 we obtain expressions which may be 
solved for. the scale factor b and the required standard deviation of 
the (y ± ) 

a 2 = In (y/X 2 ) B8 

and 

exp(-b) = Xexp(-o 2 /2) B9 

where y is the second moment of the (X^) about zero. 

Program: The log-normal generator makes use of the normal random number 

generator included in the IBM Scientific Subroutine Package for the 360 
computer. This routine (GAUSS) generates normal random numbers with any • 
required mean and standard deviation. Coding for the program is shown 
in the accompanying listing. The argument list is as follows: 

AVE - The required mean. 

VAR - The required standard deviation. 
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X - A vector of log-normal random numbers returned by the 

subroutine. Y is dimensioned by the calling program. 

^ - The number of random numbers to be generated. 

IX - A "seed" required by GAUSS. IX must be a 5 digit odd integer. 

Statements 3 to 6 compute the required standard deviation for the 
Gaus sian-random numbers and the proper scaling factor. Statements 7 to 

9 call GAUSS compute a log-normal random number from equation B5. 

Fortran List for Log-Normal Generator 

1 SUBROUTINE LOGN (AVE , VAR, Y, N, IX) 

2 DIMENSION Y(l) 

3 VAR- = VAR**2 + AVE **2 

4 SIG = AL0G(VAR/AVE**2) 

5 Z BAR = EXP (SIG/2. 0) 

6 SIG = SQRT(SIG) 

7 DO 1 I' = 1, N 

8 CALL GAUSS (IX, SIG, 0.0, X) 

9 1 V(I) = (AVE/Z BAR) *EXP (X) 

10 RETURN 
END 


11 
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